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Abstract. We report some new interesting features of the dynamics of a string 
axion field (i.e., a (pseudo-)scalar held with tiny mass with sine-Gordon-type self¬ 
interaction) around a rotating black hole in three respects. First, we revisit the 
calculation of the growth rate of superradiant instability, and show that in some cases, 
overtone modes have larger growth rates than the fundamental mode with the same 
angular quantum numbers when the black hole is rapidly rotating. Next, we study 
the dynamical evolution of the scalar held caused by the nonlinear self-interaction, 
taking attention to the dependence of the dynamical phenomena on the axion mass 
and the modes. The cases in which two superradiantly unstable modes are excited 
simultaneously are also studied. Finally, we report on our preliminary simulations 
for gravitational wave emission from the dynamical axion cloud in the Schwarzschild 
background approximation. Our result suggests that fairly strong gravitational wave 
burst is emitted during the bosenova, which could be detected by the ground-based 
detectors if it happens in Our Galaxy or nearby galaxies. 
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1 . Introduction 

The axiverse scenario provides us with an interesting possibility to hnd evidence for 
string theories through observation in the cosmological and astrophysical contexts BE]. 
String theories (M theory) are formulated in 10 (11) spacetime dimensions, and the extra 
dimensions should be compactihed to realize our four spacetime dimensions. Many 
moduli appear as a result and some of them are expected to behave as (pseudo-)scalar 
helds with tiny masses from four-dimensional point of view. Although such exotic helds 
are in hidden sector, we may detect signals from such fundamental helds through weak 
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coupling to standard model matter/fields or through coupling to gravity (see [3] for an 
overview). 

If such a scalar field with tiny mass exists in nature, what happens around a black 
hole? If the black hole is rotating, there is an ergoregion around the event horizon where 
the time-translation Killing vector dt becomes spacelike. Because of this property, the 
spatial density of the Killing energy with respect to dt of the scalar field can become 
negative in the ergoregion, although the local proper energy density with respect to a 
time-like vector is always positive if the dominant energy condition is satisfied. The 
Killing energy density becomes negative when the superradiant condition, (3.18) in 
section is satisfied. Hence, for a quasibound state mode satisfying the superradiant 
condition, the Killing energy flux towards the horizon becomes negative. Because the 
total Killing energy is conserved and represents the total energy of the system measured 
by an observer at spatial infinity, this implies that the field outside the ergoregion grows 
exponentially. This is called the “superradiant instability”, and it is one of the methods 
to extract energy from a rotating black hole analogous to the famous Penrose process. 
The superradiant instability effectively happens if the Compton wavelength of the axion 
mass n has the order of the gravitational radius of a black hole, and its typical time scale 
is about one minute for a solar-mass black hole, and about one year for an intermediate 
scale massive black hole such as the Sgr A* at the centre of Our Galaxy. The superradiant 
instability was studied also before the appearance of the axiverse scenario by theoretical 
motivations HEIEIEIIHIE] (see also uniiiimg for recent works and [13] for a review). 
The axiverse scenario has provoked renewed interests in this topic because it suggests 
that the superradiant instability may happen in our real universe, and furthermore, 
signals from an axion field may be observed by gravitational wave detectors. 

If we take account of the possibility that astrophysical black holes wear clouds of a 
string axion field, there are many issues to be studied. Up to now, some work has been 
done including simulations of scalar field as a test field [I111I5] , properties of gravitational 
wave emissions [161 imiiH], simulations including gravitational backreactions [T9|, the 
evolution of black hole parameters [THl I2U] . and so on. Our primary interest in the 
present paper is the effect of the nonlinear self-interaction of the axionic scalar field 
and its impact on the gravitational wave emission. Here, the nonlinear self-interaction 
means that the potential of the axion field becomes periodic, typically described as 

u = P/f |1 - cos(<l.//.)]. (1.1) 

This type of the potential arises by the nonperturbative effect for the QCD axion, and a 
similar mechanism is expected to work also for a string axion. In our previous paper [15] , 


we performed numerical simulations of a scalar field with the potential (1.1), i.e., the 
sine-Gordon field, as a test field in the Kerr background spacetime for a specific choice of 
system parameters. There, we numerically demonstrated that the self-interaction causes 
gradual concentration of the axion field configuration with its superradiant growth and 
eventually leads to a dynamical collapse of the configuration associated with an infall 
of positive energy towards the black hole and an outflow towards the far region. We 
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call this violent phenomenon the “bosenova” in analogy with the bosenova explosion 
observed in experiments of condensed matter physics |2l]. 


Scalar field amplitude Tbn ~ 500M 



Figure 1. Schematic picture for time evolution of the scalar cloud energy and 
emission of gravitational waves. 

Figure summarizes the overall picture of phenomena in the time evolution of 
an axion cloud with the self-interaction suggested from our previous results mm- 
In the first stage, the energy exponentially grows by the superradiant instability. At 
some energy for which the field amplitude becomes := $//a ~ 1, the nonlinear self¬ 
interaction becomes important and causes a sudden collapse of the cloud followed by 
the infall of the positive energy to the black hole and an outflow towards the far region. 
This bosenova event continues for a period ~ lOOOM, and then, the system settles to the 
superradiant phase again. Bosenova typically happens when the axion cloud extracts a 
tiny portion of the black hole rotation energy. After that, bosenova and superradiant 
growth happen alternatively. During each bosenova event, the energy extraction from 
the black hole effectively stops, and instead, a part of the axion cloud energy extracted 
from the black hole in the precedent superradiant phase is put back to the black hole. 
Thus, some part of the rotation energy extracted from the black hole is recycled. In the 
meantine, bosenova is expected to generate burst-type gravitational waves. Further, 
continuous gravitational waves are generated by oscillation of the axion field. Thus, 
the long-time evolution of the black hole is determined by the energy loss by these 
gravitational wave emissions, if the interaction with surrounding plasma and magnetic 
fields are neglected. 

In the present paper, we explore these processes more deeply. In particular, we 
focus our attention to the following three problems. First, we revisit the calculation of 
the superradiant growth rates of the quasibound states of a massive Klein-Gordon field 
(without nonlinear self-interaction). Although several calculations have been performed 
for the fundamental modes (i.e., the modes with zero radial quantum number), to our 
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knowledge no report has been made for overtone modes (i.e., modes with nonzero radial 
quantum number). Using the Leaver’s method [22], we show that there are some cases 
where an overtone mode has the largest growth rate, and thus, consists of the primary 
component of the axion cloud, depending on the system parameters and the angular 
quantum numbers i and m. In order to prove that this result is not produced by 
a mode confusion or an error, we check our results by highly accurate time-domain 
simulations. 

Next, we extend our simulations on bosenova in our previous paper na, which 
dealt with the case with a single unstable mode of a specihc type initially, to a wider 
class of initial conditions and to more complicated situations. To be specihc, we perform 
simulations for several axion mass values and mode parameters different from that in 
[TS|. For an axion cloud in the £ = m = 2 mode, our results indicate the possibility 
that the bosenova collapse does not happen in contrast to the £ = m = 1 case. Instead, 
a steady outhow of the axion held will be formed at the end point of the superradiant 
instability. We also study the ehects of nonlinearity in the case that two unstable modes 
are excited. 

Finally, we address the problem of gravitational wave emissions from the axion 
cloud. In our previous papers [16], we calculated the amplitudes of continuous 
gravitational waves emitted from this system ignoring the nonlinear self-interaction (i.e., 
from quasi-bound states of massive Klein-Gordon held) and discussed their observational 
consequences. The next step is to calculate the emission of gravitational waves from a 
dynamical axion cloud in the presence of nonlinear self-interaction. Because we already 
have the time-domain data of the dynamical evolution of the scalar held, the main task is 
to solve the perturbation equations for the gravitational held with the prescribed source 
term in the Kerr background. We attack this problem by solving the Teukolsky master 
equation in the time domain. At this moment of time, we have hnished developing 
a code to simulate gravitational waves generated by a scalar held in a Schwarzschild 
spacetime, and the extension to the Kerr case is still under progress. In the present 
paper, we report the results of test simulations for gravitational wave emissions from 
quasibound states of a massive scalar held with self-interaction around a Schwarzschild 
black hole. Although this is an artihcial setup in the sense that this system does not 
cause the superradiant instability and the initial conhguration is prepared by hand, 
we can obtain a lot of implications for gravitational waves emitted during the realistic 
bosenova. 

The present paper is organized as follows. In the next section, we briehy overview 
the basics of axions including both QCD and string axions from the ehective held theory 
point of view, and thereby explain the meaning of the basic axion parameters and the 
origin of the axion potential. Then, in section after briehy overviewing the basics 
of superradiant instability and the method to evaluate the instability growth rate, we 
show our numerical results on the growth rates of the axion cloud by superradiant 
instability taking attention to overtone modes. In section we study the ehects of the 
nonlinear self-interaction on the basis of the results of our simulations. In section we 
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explain our approach for calculating gravitational waves emitted from the axion cloud 
around a Kerr black hole, and show preliminary results of our simulations in the case 
of the Schwarzschild black hole. Section is devoted to a summary and discussions. In 
Appendix A a supplementary explanation on the pseudospectral method for calculating 


scalar held behaviour is given. In [Appendix B we give the supplementary informations 
on the gravitational wave calculation. 


2. Axions and Axiverse 

In this section, we briehy overview basic features of axion on the basis of model- 
independent arguments in order to explain the meaning of the axion parameters and 
the origin of self-interactions. 


2.1. General effective action for axions 


In this paper, we dehne an axion as a Goldstone held/boson associated with a 
spontaneous breaking of a chiral shift or U(l) symmetry. Let us denote such a chiral 
transformation with the angle parameter A by symbolically. Then, when this 

symmetry is spontaneously broken, the associated Goldstone held 0 transforms as [23] 

4 > ^ 4 > + A/a, (2.1) 

where /a is a constant called the axion decay constant in analogy with the pion decay 
constant for the pion helds as the pseudo-Goldstone helds for the spontaneous 
breaking of the approximate chiral SU(2) symmetry. The axion decay constant is 
normalized by the condition that / has the canonical kinetic term —(e/2)(V(/)^ 
(c = y/-^) and that A has the period 27r. 

A symmetry is chiral when it transforms helds with distinction of left and right 
components asymmetrically. In four spacetime dimensions, spinor helds are the only 
helds with chirality. Such a chiral transformation can be generally written 

m ^ ( 2 . 2 ) 


where t is a hermitian matrix characterizing the transformation. In higher dimensions, 
bosonic helds can have a chirality. An important example is the 5-form hux in the type 
IIB supergravity in 10 spacetime dimensions |2n [23] . 

In four spacetime dimensions, this invariance determines the coupling of the axion 
held to the spinor helds uniquely in the leading order as follows. First, because the 
action of the system is invariant under the chiral transformation, the transformation 
(2.2) of the spinor helds with hxed / produces the same change in the action as the 
transformation /—)■(/ — A/a with hxed T when A is a constant parameter. Hence, if we 
redehne T by (2.2) with A replaced by the held A = ///a, then 0 disappears from the 
terms that do not contain a derivative of T. However, the kinetic term of the spinor 
helds produces a new term as 


iTy^F.T I-)- -iTy^F.T - iTy^ 


/a 


> 1 '. 


(2.3) 
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From this, it follows that after the redehnition (2.2) of the spinor held with A = 0//a, 
the (/)-dependent part of the Lagrangian is given by 


e '£’^.0 = 


(2.4) 


if the Lagrangian does not contain a derivative term of 'L of higher dimensions. In this 
new representation, the chiral transformation does not act on the spinor held tl' and 
just shift the axion held. 

The action (2.4) is not correct in quantum theory. It is because due to the chiral 
anomaly, the ehective Lagrangian transforms as [26l EZl ESjf] 

e-'sy„ = A^i (2.5) 


a,6 


a,6 


where we have used the notation in which the coupling of the gauge held A°' = A^dx^ 
to the fermion reads 


= (a„ - iAp^) <I>, 


(2.6) 


and the gauge kinetic term can be written 


-1 




F“ = dAl“ - j/fcAl'' A A‘ 


(2.7) 


Here, it is understood that the index a runs over all the components of all gauge helds. 
So, the gauge coupling constants ga whose index belongs to the same gauge group should 
be identical. 

When we generalize the transformation to a local transformation with A = A(x), 
because of the shift symmetry of the original Lagrangian, the ehective Lagrangian 
changes as 


e-^S^eS = J^d^X + A^, 


( 2 . 8 ) 


where is the current corresponding to the chiral transformation In the path- 

integral expression for the partition function Z of the system, this change is produced 
just by the variable redehnition corresponding to the local chiral transformation [281 [2^ : 


Z 


j [dTdT ■ ■ -je'^ 


j [dT'dT' ■ ■ -je'^' 


y'[dTdT---]h(^+^^). (2.9) 


Hence, its expectation value has to vanish: 


{v„jn = (2.10) 

This modihcation of the conservation law for the chiral current can be made 
consistent with the held equations by adding the term (^//a)^ to the tree level 
Lagrangian (2.4). By this modihcation, the additional change of the ehective action 
due to the chiral anomaly is reproduced in the tree level. Thus, generalizing to the 


f In Tr{ttatb) in this expression, the left and right chirality components of the spinor fields are treated 
as independent field components. 
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multiple axion case and including the gravitational Chern-Simons (CS) term from the 
chiral anomaly, we obtain the following general effective action for axions: 

al3 a 

+ Z I (E . (2.11) 

Cf \ a,b / 

where ta is the matrix defining the chiral transformation 0 ^ is the corresponding 

axion field and = TT{tatatb)- In the multi axion case, the kinetic term of the axion 
fields may not be diagonalized if depends on other fields. In that case, the definition 
of the axion decay constant fa become ambiguous. 

Here, note that if the chiral transformations commute with the gauge 

transformations, vanishes unless both a and b belong to Abelian factors or to the 
same non-Abelian factor of the full gauge group because TT{U~HatatbU) = Trifatatb) 
for any gauge transformation U. Hence, for axions, = ^fSab except for the indices 
a and b belonging to U(l) factors of the full gauge group. If we consider a pseudo- 
Goldstone boson associated with an approximate symmetry, this rule does not apply 
and a mixing between a non-Abelian gauge field and an Abelian or non-Abelian gauge 
field can appear in the CS terms. The most famous example of such mixing is the CS 
term proportional to ti^F ■ F giving the main decay channel vr^ —)■ 27 , where F is the 
EM field. This term arises because vr^ is a pseudo-Coldstone boson associated with the 
approximate chiral SU(2) symmetry in the {u,d) quark sector. In this chiral dynamics 
approach, 7 r° corresponds to the same transformation matrix r 3 (=the Pauli matrix CT 3 ) 
as that for the third component of the weak SU(2) gauge field, and ^ = Tr(r 3 Fr 3 ) = 1. 
This produces the CS coupling of 7 r° to the EM field and the Z boson field. 

2.2. Potential and self-interactions of axion 

As we saw above, the Lagrangian of axions taking account of perturbative quantum 
corrections can be determined uniquely by the structure of the associated chiral 
transformations, or equivalently by the chiral currents ^ 5 ( 0 ), except for the kinetic metric 
and axion decay constants fa. This kinetic metric and /„ are independent of the 
axion fields but may depend on other fields and sensitive to the kinetic structure of the 
boson sector of the original theory. 

In this perturbative level, the theory preserves the shift symmetries exactly. In 
particular, there appears no mass term or potential for axions. When non-perturbative 
effects of non-Abelian gauge fields in strong coupling regions or stringy instantons 
are taken into account, however, the system acquires a non-trivial potential and as 
a consequence, the shift symmetries are broken to discrete symmetries. Although the 
exact calculation of these non-perturbative effects are out of our reach, they are usually 
estimated by the two methods [23]: the instanton approach and the effective Lagrangian 
approach for chiral symmetry breaking. 
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We first discuss non-perturbative effects of gauge field instantons. When the 
axion fields are uniform, the CS action for a G = SU(?7.) gauge field A = AHa with 
Tiitath) = Sab/2, 


Scs = epi-, pi= / ^Tr(FAF), 


( 2 . 12 ) 


is a topological quantity proportional to the 1st Pontrjagin number pi of the gauge 
bundle, where 6 is expressed in terms of some constant 6q and the axion fields as 

e 


0 = 0o + ^ Tr(4ta4) = 


(2.13) 


In the Lorentzian spacetime, there exists no classical solution for which pi is bounded, 
while in the Euclidean space obtained by the Wick rotation, there may exist a classical 
solution with finite pi. Such a solution is called an instanton and give a non-trial 
contribution in the Euclidian path integral formulation. 

Topologically, instanton solutions in four dimensions are classihed by 7r3(G). 
Because 7r3(U(l)) = 0, there exists no instanton for an Abelian field, hence the CS term 
induces no potential. In contrast, because 7r3(SU(n)) = Z(n > 2), all non-Abelian SU(n) 
gauge field has instanton solutions topologically classified by the instanton number 
u G Z. Among these solutions, those with self-dual or anti-self-dual field strength have 
the least action for a given instanton number, hence give the dominant contribution in 
the path integral. Explicit solutions are known for some cases. The most famous one 
is the self-dual BPST solution with pi = I for the SU(2) gauge field [20], which is also 
an instanton solution for the SU(n) gauge field with u > 2 by the standard inclusion 
SU(2) C SU(u). Because the BPST solution has the moduli freedom corresponding to 
the spacetime translation and the dilation, its contribution to the partition function Z 
can be expressed as 


Zi = A^ / d^x / dR e 


—S'b 


(2.14) 


where A is some mass scale and R is the size of the instanton. For the self-dual 
solution, the Euclidian instanton action S'b is determined by the coupling constant 
and the instanton number as 


Ae= [ d^x\TT{*F AF) = ±\ [ d^xTr{F AF) = ^\pi\. 
J 9 9 J 9 


(2.15) 


If we assume that a general instanton solution with pi = n can be approximated 
by the linear superpositions of p BPST solutions and of q anti-BPST solutions with 
n = p — q [the dilute gas approximation), the total contribution of the instantons to Z 
can be written 


^ n\ q\ 
P,<1>0 


A*' / d*x / + 


.(2.16) 


Hence, the instantons produces the non-perturbative potential 

V = —A^cos6' -|- const. = —A^cos j -|- 'y^C,°‘4>a/fa 1 + const., (2.17) 
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where 

= 2A® J dR (2.18) 

From this formal expression, we see that instantons with large size R produces 
large contribution if the gauge field becomes strongly coupled in the IR region, but we 
cannot calculate this strong coupling effect properly. In the QCD axion case, A is of the 
order of the pion mass as we see below. In the case in which there exists more than one 
gauge fields with strong coupling regions, the axion potential is given by the sum of the 
contributions from all such gauge fields. Thus, generically, the potential is given by 


R = - ^ A^ cos j + COUSt. (2.19) 

Here, OAflA are phases determined by the CP phase of the fermion mass matrix and the 
^-angle of the vacuum of each non-Abelian gauge field. 

Note that if the number of independent axions is smaller than the number Ns of 
strongly coupled gauge fields, some of Oa^A may not be set to zero by the constant shift 
of the axion fields. In this case, CP invariance is not recovered in some gauge sector 
by the Peccei-Quinn type mechanism. In particular, a strongly coupled gauge field in 
the hidden sector may interfere the QCD sector to cause a CP violation larger than the 
experimental upper limit. In contrast, in the case > Ng, 6 a,o can be always set to 
zero by shifts of the axion fields as far as the effects of instantons of four-dimensional 
gauge fields are taken into account. However, these residual shift symmetries may be 
broken by stringy effects as discussed below. 

The masses of the axions can be calculated by expanding V with dA,o = 0 with 
respect to (pa and diagonalizing the quadratic terms 


R2 = 



( 2 . 20 ) 


From this expression, we see that the axion mass is roughly given by 

N 

— ( 2 . 21 ) 

/a 

where A is the largest strong coupling energy scale of the gauge fields coupled to the 
axion and /a is the axion decay constant. 

Another important consequence of this formula is that if N^, > Ng, N^, — Ng axions 
should be massless. This problem is serious for the axiverse scenario because Ng cannot 
be become so large. We can easily show that the kinematic mixing of the axions does 
not change the situation. This problem can be avoided if we take into account string 
instanton effects [311132] 

The second method to estimate the non-perturbative effects is to use the effective 
field theory for composite mesons corresponding to the pseudo-Goldstone bosons 
associated with a spontaneous breaking of non-anomalous approximate chiral symmetry. 
For example, let ns consider the chiral SU(2) symmetry of the QCD system with two 
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quarks u and d. First, let us move the 6 phase (2.13) for the color SU(3) gauge field 
to the quark mass matrix sector by the chiral transformation [u, d) —)■ 
with Uu + Ud = Then, the quark mass terms changes as 

imuUU + iniddd —)■ d, (2.22) 


and the kinetic term of the quarks produces the additional axion quark coupling 


df,9{yuU'y^j5U + ydd-f^'j^d). 


(2.23) 


Now, in the strong coupling regime, the products of quark fields acquire non-vanishing 
vacuum expectation values as 


— i (uu) = —i (dd) = UcCos (27r°//,r) , (2.24a) 

- i (m 75 m) = i {d'y5d) = -iUc sin (27r°//^) , (2.246) 

(h7^75M) = - (d7^75c() = (2.24 c) 

where 7r° is the field of the neutral pion and /t^ is the pion decay constant. Then, the 
quark mass term produces the potential 

VaTT = -VcTUu cos {yu9 - 271°/f^) - VcTUd COS {yd9 + 277°/ f^) , (2.25) 


and the derivative coupling between the axions and quarks gives rise to the axion-vr® 
kinetic mixing 


1 

2 



4) + ^^^d{yu - yd) 




(2.26) 


where we have set taU = z°u and tad = z°[d. 

In the case of a single axion, which is supposed to be the QCD axion, the value 
of {yu,yd) is uniquely determined by the requirements that yi + y 2 = —1 and that this 
axion-pion mixing disappears. Then, from 14^, we obtain the standard formula for the 
axion mass and the pion mass in terms of Vc- 


m. 


~ AVr 


m,. 


rud 


ml ~ u, 




f2 ’ 

J TV 

rriumd 


/2 mu + md 


a, 

2/a 


mumd 


{mu + md)2 




(2.27a) 

(2.276) 


where we have assumed the invisible axion condition fa,^ fw- 

In contrast, if two or more axions couple to the color SU(3) gauge fields or the 
quarks u and d, we cannot eliminate the vr^-axion mixing by the choice of yu and yd- 
This mixing, however, can be eliminated by just shifting the axion fields other than the 
QCD axion by some constant multiple of 77° field. Thus, the axion-pion mixing does 
not provide additional mass to axions other than the QCD axion, and just yield a tiny 
correction to the 77° mass. Actually, we can show that the kinetic mixing between axions 
and mesons does not change the number of massless axions [33] . 
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Figure 2. Various cosmophysical phenomena that may be caused by string axions 
and constraints by future experiments in the plane of the axion parameters (rua,/a)- 
The CMB experiments PIXIE [34] /PRISM [35] and the axion helio scope experiment 
lAXO [36] utilize the CS coupling of an axion to the EM field, g^jE ■ B. Hence, they 
do not directly give constraints on In this figure, we have translated the expected 
sensitivity on to the constraint on /a by the relation g^,^ « a/{'!TfA- 


2.3. Axiverse 

In string theory, the compactihcation of extra dimensions can produce a large number 
of axions if the internal space has a rich topological structure. This is because the 
generalized gauge symmetry for form helds contained in the original ten or eleven 
dimensional supergravity turn to a shift symmetry of zero modes. As an example, let 
us take up the RR 2-form held C 2 = {l/2)CMNdx^dx^ in the type IIB supergravity in 
ten-dimensional spacetime. When it is compactihed to the product of a four-dimensional 
spacetime A 4 and a six-dimensional Calabi-Yau manifold Ig, the zero mode of C 2 can 
be expanded in terms of harmonic 2 -forms X 2 ^6 up to the gauge transformation 
C *2 —t C *2 T dcT\ as 

C2^oj2{x)^^(t)^{x)x2. (2.28) 

a 

where uj 2 {x) and (l)a{x) are a 2-form and functions on X 4 , respectively. Because 
a closed form is locally exact, the theory is also invariant under a transformation 
C 2 —)■ C 2 + /d 2 , where 132 is an arbitrary closed form. This implies that the theory 
is invariant under the shift of the helds, t 0 a (a^) -|- Ca with Ca being constants. 

Further, we can conhrm that the term of the form ^ appears in the 

ehective four-dimensional theory after compactihcation, where 'ip/s are spinor helds on 
X 4 corresponding to the zero modes of the gravitino held in the internal directions 
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and the dilatino field on Yq. This implies that the current of this shift symmetry 
•^ 5 ( 0 ) ~ '^j k is chfial, and 0 q’s become axion fields according to the general 

arguments in the previous subsections. Thus, we obtain 62 (^ 6 ) axions from C 2 , where 
62 is the second Betti number. 

In the original type IIB theory, there exists no matter gauge fields, and no CS term 
appears. However, if we introduce D-branes, the CS D-brane action 1371 



(2.29) 


produces the CS interaction of the axions to gauge fields and the gravitational field on 
the branes after compactification. This is natural because this CS action is obtained 
from the anomaly free condition. 

All string axions should be described by the effective theory in four dimensions 
explained in the previous subsections. In particular, their roles in cosmophysical 
phenomena and experiments are essentially determined by its mass rua and axion decay 
constant fa{see figure § 0 . We can probe the compactification structure and string 
theory behind it mainly through the distribution of the values of these quantities. 

From the perspective of searching phenomena that may be caused by axions, it is 
important to know plausible predictions of theory on the parameter distribution. As we 
mentioned in the previous subsection, we can give non-perturbative mass by instanton 
effects of low energy gauge fields only to W axions where Ng is the number of gauge 
fields with strong coupling regimes including those in the hidden sector. Hence, in the 
axiverse perspective, most of string axions get mass by string or brane instanton effects. 
In this case, in the mass formula ma ~ A^//a, /a is generally given by /a ~ M^i/S e, and 
A^ ~ where S'e is the instanton Euclidean action [321 EH. Svrcek and Witten 

estimated S'e by considering the interference of string axions to the QCD theta phase 
and obtained the constraint S'e > 200 , i.e. /a < lO^^GeV. On the basis of this result, 
the axiverse paper [ 1 ] argued that mg. is distributed uniformly with respect to logma 
because nig oc 

As we have shown in figure various experiments that may detect axions or provide 
constraints on the axion parameters are planned. These experiments, however, require a 
sufficient abundance of axions or rather strong coupling of axions and photons, hence a 
rather small value for the axion decay constant. In contrast, the superradiant instability 
of an axion field around a rotating black hole that we discuss in the present paper does 
not require any abundance because the seed can be zero point oscillations of the axion 
field. Further, the phenomena caused by the instability becomes more violent and 
energetic for a larger fg as we will see. Hence, this phenomenon enables us to probe the 
parameter region complimentary to other experiments. 
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3. Superradiant instability 

In this section, we study the growth rate of the superradiant instability for a massive 
Klein-Gordon held without nonlinear self-interaction. We focus on overtone modes. 


3.1. Formulation 


We consider a Kerr black hole whose metric is given by 

, n A — sin^ 9 . n AMar sin^ 9 , , 
ds^ =-^-df^--dfd(^ 


S S 

''(r^ -|- a?Y — Aa^ sin^ 9 


sin^ ^-dr^ + Sd0^ 


with 


'F = cos^ 6*, A = — 2Mr, 


(3.1) 


(3.2) 


in the Boyer-Lindquist coordinates. Here, M is the black hole mass, and the angular 
momentum is given by J = Ma. We often use the nondimensional rotation parameter 
a* := a/M. The event horizon is located at r = r+ := M -|- — a?, and it rotates 

with the angular velocity Hh = ®/(^+ + The important property of the Kerr 
black hole is the existence of the ergoregion where the time coordinate basis = {dtY 
becomes spacelike, i.e. = gu > 0. Because of this property, the Killing energy of a 
particle/boson held with respect to dt can be negative in the ergoregion. 

In this section, we study the behaviour of massive Klein-Gordon held, whose 
behaviour is governed by the massive Klein-Gordon equation. 


- /i^G' = 0 

with U = The separation of variables is possible in this case as 

$ = 2Re [e-^‘^^R{r)S{9)F^^] . 

The equations for S{9) and R{r) become 


1 d 
sin 9 d9 


. AS 


cos^ 9 


m 


sin^0 


+ A. 


£m 


s = o, 


dr 


dR\ ) [(r^ -I- a^)uj — am]" 


^ + 


dr / 


A 


A 


£m 


R = 0, 


(3.3) 

(3.4) 


(3.5) 

(3.6) 


where 


; 2 2 2 
k = jj, —CO , 

Aim = Aim + a^ca^ — 2amu. 


We normalize each angular mode function S as 


d^sin^l^p 


1 


(3.7a) 

(3.76) 

(3.8) 
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In order to study the behaviour of R near the horizon, it is useful to introduce the 
tortoise coordinate r* by 

dr. = (d^dr. (3.9) 

In this coordinate, the horizon is located at r* = —oo, and from the regularity of 
<1> around the future horizon in the future Finkelstein coordinates, we hnd that the 
behaviour of the radial function R{r^) becomes 


for ingoing waves. 

Our problem is as follows. We seek solutions that fall off at inhnity and satisfy the 
ingoing boundary condition (3.10) near the horizon. Since the two boundary conditions 
are imposed, the problem becomes an eigenvalue problem. There is a close analogy 
between our case and the hydrogen atom in the quantum mechanics. In the hydrogen 
atom, the regularity at the origin and the falling-off condition at inhnity are imposed, 
and as a result, the energy levels are discretized. Similarly, in our problem, the angular 
frequency u is discretized. Since there is a hux towards the horizon, the eigenfrequency 
becomes a complex number in general. 


a; = cjR -|- icui. (3-11) 

If the imaginary part Ui is negative, the held decays in time, and the system is stable. 
In contrast, if ui is positive, the held grows exponentially, and the system is unstable. 
This instability is called “superradiant instability”. 

In order to understand the mechanism of the superradiant instability, it is helpful 
to introduce the Killing energy E{t). The energy-momentum tensor of a scalar held is 

(Vp$V^<l> + 2Um , (3.12) 

and the conserved energy current is introduced by 


pM = (3.13) 

After some algebra, the energy contained in the part r* > of the t = const, 
hypersurface is found to be 


£(i) = 


h=const,r(‘"^ 


P'^dS, 


1 

2 


„(in) 


{r^ + a^)^% + 


EA 

-f 


+ 2C/(4) 


(3.14) 


dr*df2. 


with dfl = sin0d6'd0 and $ = d<h/df. On the other hand, the energy hux towards the 
horizon with respect to the time coordinate observed at r* = is 


Mt) = 


I 


<I),^>(r2 + a2)dO. 


(3.15) 


For these formulas, the energy conservation holds: 


E Ee — 0. 


(3.16) 
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-— K, Fe ^ 2u{u: — 'mVti{){r\ + a^). 


(3.17) 


Now, we direct our attention to the near-horizon region and assume that the held 
behaves as (3.10). The Killing energy density with respect to the tortoise coordinate 
di?/dr* and the hux towards the horizon, Fe, are 
dE 
dr* 

Therefore, if the superradiant condition 

u < mfln (3.18) 

is satished, negative energy distributes in the neighborhood of the horizon and it falls 
into the black hole. Then, because of the energy conservation (3.16), the energy outside 
of the horizon increases. This explains the origin of the superradiant instability, i.e. an 
unstable mode with a positive ui. 

For a hxed a*, the system can be specihed by a nondimensional parameter M/r if we 
measure lengths in the unit of M. For M/x ^ 1 and M/x 3> 1, the eigenvalue problem 
can be approximately solved by the matched asymptotic expansion method [1] and by 
the WKB method [5], respectively. In the matched asymptotic expansion method for 
Mfi -C 1, the equation is solved in the near-horizon region and in the far region, and 
they are matched to each other in the matching region. The solution in the far region is 
the same as that of the hydrogen atom in quantum mechanics, except that is replaced 
by Mjji. Therefore, a quasibound state of a massive scalar held can be regarded as a huge 
gravitational atom, and the solutions are labeled by the angular quantum numbers ^ 
and m, and the radial quantum number n-^ = 0,1, 2,... that characterizes the oscillation 
in the radial direction (or equivalently, the principal quantum number Up = £ -|- 1 4- n^). 
In this paper, we call the solution with Ur = 0 the fundamental mode, and the solutions 
with Ur > 1 the overtone modes. 

In the parameter region Mfi ~ 1, this eigenvalue problem has to be solved 
numerically, and some work has been done on this problem illSlE]. Here, we follow 
the method by Dolan |9j, where the radial function R{r) is assnmed to have the form 
of the inhnite series. 


R{r) = (r — r+) “^(r — r_) 


icr+X-^a-kl' 


n=0 


with 


2Mr i 


a = 


(OJ 


mQ 


H 


and X = 


M{2uj^ - fR) 


(3.19) 


(3.20) 


r+ — r_ ' ' " k 

Here, Re(fc) > 0 is required in order to make the radial function fall off at the distant 
region, r ^ M. This expression correctly reproduces the held behaviour near the horizon 
and at inhnity. Substituting the expression (3.19) into the radial equation (3.6), we have 
a three-term recurrence relation for a „: 


_ /3o 

Oi —-Oo, 

T l3n(^n T 1 0) 


(3.21a) 

(3.216) 
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where 


= (n + 1) [n + 1 — 2ia], 
I3n = — ‘2n^ + 2n 


M 


— 1 + 2icr — 2bk -— 2a;^) + 

k 


2M, , p, 

n + 1 — 2icT-( 20 ;^ — /i ) 

k 


In = {n- 1) 

with b = y/M"^ — a? and 

= a^k"^ - 2M{M + b){n^ - 2a;2) + - (1 - 2iMuj) 

M 


+ 7n 


(3.22o) 

(3.226) 

(3.22c) 


1 + -iam — 2 M‘^oj) 
b 


-2A; - (/i^ - 2n;2) 

bk 


[6(1 — 2iMa;) + i(am — 2M^a;)] — (3.23a) 


7i = 


- 2a;2)2 - 2uj^) 




bk 


[—am + 2M'^oo + 2i6(l — iMoj)\ 


+ (1 - 2iMa;) 


1 + - (am — 2M'^u) 


(3.236) 


Although the appearance of these formulas differs from that of [S], it can be checked 


that they exactly agree. From the three-term recurrence relations (3.21a) and (3.216), 


we obtain the following equation for u in the form of the continued fraction: 

0 = (3.24) 

Pi— P2— P 3 — 

which can be solved numerically using the Newton-Raphson method. In the numerical 
computation, we took account of the first 10^ terms. One subtlety in this method is 
that we have to calculate the angular eigenvalue A^m simultaneously. We avoided this 
problem by using the approximate formula for small ak derived by Seidel [3H] (see also 
|39j). which is the series expansion up to sixth order. Seidel’s approximate formula 
gives fairly accurate value in the problem below, and the error from the truncation of 
the series expansion is of O(10“^a;i). 

In our previous paper [3], we developed a numerical code to solve this problem, 
and some results are shown in mm- In the present paper, we report numerical results 
for overtone modes obtained by the same numerical code. As far as we know, all works 
up to now took attention just to the fundamental mode = 0, and no results for the 
overtone modes have been reported. 


3.2. Numerical results 

The overtone modes are obtained by just changing the initial guess in the Newton- 
Raphson method. To be more specific, we first choose a value of M/i for which 
the imaginary part of the quasibound state solution is small. Then, we solve the 
equation for u with various initial guess for a;//i, which is taken to be a real number, 
by gradually increasing from a small value to unity. Typically, a quasibound state 
solution with a larger number of n,. has a larger real part of the eigenfrequency cur. 
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After several solutions are obtained, we determine the radial quantum number by 
drawing the wave function explicitly. The sequence with a fixed n^- can be obtained 
by repeating the process adopting the solution for some Mp as the initial guess for a 
slightly larger/smaller value of Mfi. 

Although the Leaver’s method has the ability to determine the solution with 
arbitrary accuracy in principle, in some cases the solution becomes inaccurate by hitting 
the machine accuracy. This tends to happen as or a* is increased and Mfi is decreased. 
Such data are omitted from the figures below. 


at = 0.90 at = 0.99 




at = 0.999 



Mfl 

Figure 3. The imaginary part of the boundstate frequency, Mwi (i.e., growth rate of 
the superradiant instability), for = 0, 1, 2, 3, 4, and 5 (red, green, blue, purple, sky 
blue, and grey, respectively, online) and ^ = m = 1, 2, and 3 for the case a* = 0.90 
(upper left), 0.99 (upper right), and 0.999 (bottom). Each inset highlights the peaks 
of six curves for £ = to = 3. 

Figure shows the values of the imaginary part of the boundstate frequency, Mwi, 
for a* = 0.90 (upper left), 0.99 (upper right), and 0.999 (bottom). In each panel, the 
results for the modes with £ = m = 1,2, and 3 and n-^ = 0, ...,5 are shown. In the 
case of £ = m = 1, the growth rate of Ur = 0 is greater than that of Ur = 1, and they 
are different by a factor of a few. For a* = 0.99 and M/i = 0.40, for example, the two 
growth rates are different by a factor of 2.6. The difference between the two growth 
rates is smaller for i = m = 2: For a* = 0.99 and Mfi = 0.80, for example, the difference 
is 22%. 
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An especially interesting thing happens for i = m = 3, as shown in the inset of each 
hgnre. In this case, there is a parameter region of Mfi for which the growth rate of the 
overtone mode with n,- = 1 becomes larger than that of the fnndamental mode n,- = 0. 
This happens in the domain 0.771 < M < 0.920 for a* = 0.90 and in the domain 
0.772 < M/i < 1.366 for a* = 0.99. It is also worth mentioning that for a* = 0.99, 
there are regions where the growth rates of the = 2 and 3 modes are larger than that 
of the fnndamental mode rir = 0, althongh in these regions the most unstable mode 
is the Ur = 1 mode. When the black hole is close to an extremal one, there appears 
a domain where the higher overtone mode with = 2 has the greatest growth rate: 
For a* = 0.999, this happens for 1.353 < M/x < 1.508. On the other hand, the growth 
rate of the = 1 mode becomes largest in the domains 0.773 ^ M/i < 1.352 and 
1.509 < M/i < 1.564. In the other domains, M/i < 0.772 and 1.565 < Mp, < 1.656, the 
nr = 0 mode has the largest growth rate. 

Up to now, most discussions on the evolution of axion cloud have been done under 
the assumption that the Ur = 0 mode is most unstable. However, our result indicates 
that when the black hole is rapidly rotating, an overtone mode may grow hrst in the 
£ = m = 3 mode. Since this result is very important, we will check its validity using a 
different method below. 

a* = 0 . 99 ,1 = m = 1 a* = 0 . 99 ,1 = m = 2 




r*/M r^/M 


a* = 0 . 99 ,1 = m = 3 



r^/M 


Figure 4. Snapshots of the quasibound states of the modes t = m = 1 (upper left), 2 
(upper right), and 3 (bottom) for Mp = 0.40, 0.80, and 1.30, respectively. The black 
hole rotation parameter is a* = 0.99. Five curves in each figure correspond to = 0, 
1, 2, 3, and 4 (red, green, blue, purple, and sky blue, respectively, online). 
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Figure l^shows snapshots of the wavefunctions for £ = m = 1 (upper left), 2 (upper 
right), and 3 (bottom) in the case a* = 0.99. The values of M/i are 0.40, 0.80, and 
1.30, respectively, and the modes with radial quantum numbers = 0, ...,4 are shown. 
The field configuration for r^/M > 0 is very analogous to the wavefunction of the 
hydrogen atom: The n,- = 0 mode has one local maximum, the = 1 mode has one 
local maximum and one local minimum, the = 2 mode has two local maxima and 
one minimum, and so on. The position of the first peak slightly becomes closer to the 
horizon as is increased, and the mode with larger rir tends to distribute also at the 
distant region. 


3.3. Time evolution for i = m = 3 


Because the result for the case i = m = 3 may be counter-intuitive, we have carefully 
checked the correctness of the numerical method. First, we compared our numerical 
data with those obtained by Dolan [9] for the case Ur = 0. The two results agree very 
well. Next, we calculated the time evolution of the quasibound states for the situation 
where the growth rate of Ur = 1 is larger than that of Ur = 0, starting with the initial 
data generated by the Leaver’s method. As an example, we chose the system parameters 
as a* = 0.99 and Mjj, = 1.35. This provides ns with a good test because if the result 
of the Leaver’s method is generated by some mistake or error, the field configuration 
would include some spurious data, and then, the field would relax to a different state 
and the growth rate would change in the time evolution. 

For this purpose, we used our code based on the pseudospectral method described 
in the next section. Our code is sixth-order in the radial direction and fourth-order 
in the time direction, and the grid differencing with Ar* = 0.2M and At/Ar* = 0.05 
turned out to be sufficient to achieve required accuracy. In order to evaluate the growth 
rate in the time-domain simulation, we used the two formulas. 


Ut — 


and = 


E 


(3.25) 


2E{ty “““ 2E{t) 

Here, E{t) denotes the total energy from r* = —lOOM to r* = 300M, E := di?/df, and 
Ee denotes the flux towards the horizon at r* = — lOOM. The latter formula for a;) 
holds because the energy is quadratic in the field $ for the Klein-Gordon field, and E is 
evaluated from the numerical result of E{t). The former formula for holds because 
of the energy conservation (3.16). Although the energy conservation must hold exactly, 
they slightly differ in the numerical simulation because of the numerical error. 

Figure shows the result for the time evolution up to f = lOOOM. The left panel 
shows the energy density dE/dr^ with respect to the tortoise coordinate at time t = 0 
(thick grey lines) and lOOOM (thin black lines) for the cases = d and 1. For each of 
Ur, the energy distribution scarcely changes and we interpret this result as the evidence 
for the absence of a spurious contribution from numerical errors. Next, let us look 
at the energy amount more closely. The right panel shows the increase in the energy 
/S.E ;= E{t) — F^(0) normalized by the initial amount of the energy F^(0). The numerical 
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a, = 0.99, Mn= 1.35 a, = 0.99, M^l= 1.35 



Figure 5. Results for the time-domain simulation of the quasibound states in the 
modes £ = m = 3 and = 0 and 1. The system parameters are a* = 0.99 and 
Mfi = 1.35. Left panel: Comparison of energy density with respect to the tortoise 
coordinate r^,/M at t = 0 (thick grey lines) and t = lOOOM (thin black lines) for 
Ur = 0 and 1. The two lines approximately coincide for each Ur- Right panel: Time 
evolution of increase in the energy AE normalized by the initial total energy E{0) for 
Ur = 0 and 1. The energy increase rate for Ur = 1 is larger compared to that for nj. = 0. 


data of the energy E{t) show a clear signal of growth in time, and the growth rate for 
nr = 1 is larger than that for = 0. 


Table 1. Comparison between the growth rate Mlui of the mode i = m = 3 and 
Ur = 0 and 1 for the case M/r = 1.35 and a* = 0.99, evaluated by the three different 
methods. They are indicated by and (see text for definitions). All the 

three methods give consistent values for each of nr = 0 and 1. 


fir 




0 

3.1763 X 10-9 

3.1754 X 10-9 

3.1736 X 10-9 

1 

4.4789 X 10-9 

4.4784 X 10-9 

4.4758 X 10-9 


Table summarizes the growth rate obtained by the Leaver’s method, and 
and by the time evolution. During the time evolution, and approximately 
stay constant, with typical fluctuation of 10“^^ and 10“^^, respectively. All of these 
values are consistent for each of rir = 0 and 1. This result is the evidence for the 
correctness of the results by the Leaver’s method: For i = m = 3, an overtone mode 
can have a greater growth rate compared to the fundamental mode when the black hole 
is rapidly rotating. 

This result may seem strange at hrst. However, when the radial quantum number 
TT-r is increased, the position of the hrst peak shifts towards the horizon as sho wn i n 
hgure|^ and therefore, the tunneling effect becomes stronger. In terms of of (3.25), 
the numerator —Fe becomes larger if we compare with a hxed value of the hrst peak. 
On the other hand, since a larger amount of the energy distributes at the distant region 
for larger the denominator 2E{t) also becomes larger. Therefore, these two ehects 
compete with each other, and it is nontrivial whether uJi decreases or not when is 
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increased, especially when M/a is large. 

3.4- Summary of this section 

In this section, we revisited the calculation of the superradiant instability growth rate 
for quasi-bound states of the massive Klein-Gordon held taking account of the overtone 
modes. For i. = m = 1 and 2, the growth rate decreases by a factor when the radial 
quantum number n^- is increased by one. For i = m = 3, there appears a parameter 
region of Mfi where the overtone mode with = 1 gives the largest growth rate rather 
than the fundamental mode with = 0. If the black hole is spinning very rapidly, 
a* = 0.999, there appears a region where the mode with Hj- = 2 gives the largest growth 
rate. 

Our results suggest the importance of overtone modes in studying the evolution 
of the axion cloud around a rotating black hole, especially when the angular quantum 
numbers i and m are large. Most of the studies on the axion cloud up to now are based 
on the assumption that the most unstable mode is the fundamental mode n,. = 0. For 
example, in our previous work [12], we calculated the gravitational wave emission rate 
from the axion clouds in the modes i = m = 1,2, and 3, and we used the quasibound 
state solutions of the fundamental mode rir = 0 also for i = m = 3. In order to make the 
calculation more realistic, a reconsideration using the most unstable mode is required, 
although we postpone it as a future problem. 

Also, when i and m are one or two, there is a possibility that the overtone modes 
contribute to the axion cloud since the typical time scale is not very different from 
that of the fundamental mode. In the next section where the effect of nonlinear self 
interaction of axion cloud is studied, we also take attention to the axion cloud that 
consists of superposition of the modes Ur = 0 and 1. 


4. Axion Bosenova 


In this section, we study the non-linear evolution of a massive scalar hied due to self¬ 
interaction in the Kerr-background spacetime. We adopt the equation (3.3) but with the 
potential U replaced by (1.1). It is quite useful to normalize <F with the decay constant 

/a, 

(4.1) 


The equation becomes the sine-Gordon equation, 
V‘^(p — sin 9 ? = 0. 


(4.2) 


When (p <C 1, the held behaviour is well described by the massive Klein-Gordon 
equation, and the simple superradiant instability discussed in the previous section takes 
place. However, when ip grows to order unity, the ehect of nonlinear self-interaction 
becomes important. 
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In our previous paper na, we reported only the result for the specihc quasibound 
state with i = m = 1 and n,- = 0 for the system parameters a* = 0.99 and M/i = 0.40. In 
this section, we show the results of our simulations for a wider class of initial conditions 
and to more complicated situations including mutimode excitations. 


4-1. Numerical method 

The three-dimensional (3D) code to simulate the dynamics of a scalar held in the Kerr 
background was developed in our previous paper na. Since then we also have developed 
another code to simulate the same system based on the pseudospectral method that is 
sometimes very powerful as we have seen in the previous section (but depending on the 
situation). We explain these codes briehy one by one. 


4-1.1. 3D code The 3D code solve (4.2) by the hnite differencing method both in the 


radial and angular directions. The subtle point is that the numerical instability rapidly 
develops if we use the Boyer-Lindquist coordinates (f,r*,6*,0). Our prescription is to 
adopt (f) dehned by 


0 O^amoh 


(4.3) 


instead of 0. Here, Dzamo is the angular velocity of the zero-angular-momentum 
observers (ZAMO) staying at hxed values of r* and 6 and rotating in the (j) direction 
keeping vanishing angular momentum. With this coordinate, stable simulations become 
feasible. The interpretation is that the frame associated with the Boyer-Lindquist 
coordinates propagates superluminally in the ergoregion, and the numerical evolution 
equation cannot respect the causality. On the other hand, the ZAMO frame are 
timelike everywhere and the causality is properly included in the evolution equation. 
Since the ZAMO coordinates gradually become distorted, we pull back the coordinates 
periodically. See [12] for more details. 


4-1.2. Pseudospectral method Our new code is based on the pseudospectral approach, 
where the angular dependence of the function (p is spectrally decomposed. First, we 
decompose the coordinate 0 as 

(p= Y, sinl"*l0/(”^)(t,r,0)e™‘^, /(-»«) =/My (4.4) 

m=0,±l,±2,... 

The latter condition is imposed in order to guarantee that is a real number. As for 
the 9 direction, we use y dehned by 

y = cos 6* (4.5) 

rather than 6, and decompose the coordinate y as 

/<”)((. Y. (4.6) 

n=0,l,2,... 

The expansion of this form is practically same as the expansion by the spherical 
harmonics. Then, we solve the equation for a^\t,r^) by the hnite diherence method 
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with the sixth order in the r* direction and the fourth-order Runge-Kutta method in the 
time direction. In the numerical simulation, we introduce the cut-off parameter rimax 
and mmax, and require = 0 for n > rimax and m > rri ma x. 

In order to derive the equation for , we have to express sin (p also in the form 

(4.7) 


with 


sm(p= sinH0sM(t,r,0)e^™'^, 

m=0,dil,di2,... 

n=0,l,2,... 


g(-m) _ g{rn)* 


(4.8) 


The calculation of can be reduced to the algebraic computation using the Taylor 
series of sin(p. The details of this method and the equation for are presented i 
Appendix A 


m 


We note one subtlety in this decomposition. Due to the reality condition = 


(4.4), the distinction of the quantum numbers between -|-m and —m apparently 
loses a meaning. However, it is not the case because contains both positive 
frequency mode (cx and the negative frequency mode (cx e'‘^*) with a; > 0. 

Then, the positive and negative frequency modes are naturally reinterpreted as the 


real -|-m and —m modes, respectively. Therefore, introduced in (4.4) represents 


the superposition of the ±m modes. In explaining the simulation results below, we often 
refer to the excitation of the m = — 1 mode defined in this way. 

Here, we discuss the features of this new code. Since the convergence with respect 
to m is relatively fast, the simulation is not as heavy as the 3D code. Also, the 3D 
code requires a small time step for a stable simulation, At/Ar* = 0.05, while the 
pseudospectral code allows relatively large time step, Af/Ar* = 0.25, although a small 
time step must be adopted for a calculation that requires a high accuracy like the 
one in the previous section. In the case where nonlinear self-interaction is present, 
the calculation of causes a large number of operations, and the simulation time 
becomes much longer compared to the massive Klein-Gordon case. But in general, the 
calculation of the pseudospectral code is much faster compared to the 3D code. 

The numerical stability is a little problematic. In the pseudospectral approach, 
a stable simulation is available as long as the modes with large n remain small, and 
the results by the two codes are consistent. But when violent dynamics causes a huge 
excitation of the modes with large n, the simulation easily crashes. Although simulations 
of bosenovae could be done in some cases, we found many cases where a bosenova causes 
a crash of the simulation. However, since the simulation is faster compared to the 3D 
code, it can be used as a tool to search for the system parameters that lead to interesting 
phenomena. In this sense, having the two codes is fairly useful for our investigations. 


4 . 2 . Simulations 

Now, we present the results of several simulations including nonlinear self-interaction 
effects. The setups are summarized in Table In all simulations, we fix the black hole 
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rotation parameter to be a* = 0.99. We choose the initial condition to be (sum of) 
the quasibound state solutions in the Klein-Gordon case (in the absence of nonlinear 
interaction). Then, the system parameter is specihed by M/i and the amplitudes of the 
modes that are labeled by the angular quantum numbers {£, m) and the radial quantum 
number n^. 


Table 2. Setup of the performed simulations in this paper. 


Simulations 

Initial {£,m,nA [Amplitude] 

Mu 

Bosenova? 

(la) 

(1,1,0) [0.40] 

0.30 

No 

(lb) 

(1,1,0) [0.45] 

0.30 

Yes 

(2) 

(2,2,0) [1.00] 

0.80 

No (but see text) 

(3) 

(1,1,0) [0.70] and (2,2,0) [0.01] 

0.40 

Yes 

(4) 

(1,1,0) [0.60] and (1,1,1) [0.20] 

0.40 

Yes 


In section 4.2.1, we present the simulations (la) and (lb), where the initial condition 
is chosen to be a pure mode with £ = m = 1 and Ur = 0 for Mfi = 0.30. The same 
situation was studied in our previous paper [T3] but for Mfi = 0.40. Here, we are 
interested in the dependence on the axion mass, especially how the condition for the 
bosenova occurrence changes with the mass. In section 4.2.2 , we present the results for 
the case (2) where the initial condition is represented by a pure mode with £ = m = 2 
and Ur = 0. 

In simulations (3) and (4), we study what happens if the initial condition is given 
by superposition of two modes. In section 4.2.3, we consider the case (3) that initially 
the axion cloud is dominated by the £ = m = 1 mode but a small contribution of the 
£ = m = 2 mode is present. This situation is realistic because for the chosen mass, 
Mfi = 0.40, both the £ = m = 1 and 2 modes are unstable and the growth rate of the 
£ = m = 2 mode is much smaller than that of the £ = m = 1 mode. Here, we are 
interested in whether the £ = m = 2 mode can continue to grow or not in the situation 
where the nonlinear effect of the £ = m = 1 is strong. In section 4.2.4[ we consider the 
superposition of the fundamental and overtone modes (nj. = 0 and 1, respectively) in 
the £ = m = 1 case and study to what extent the property of the bosenova is changed. 


4 . 2 . 1 . Pure £ = m = 1 mode We begin with the cases (la) and (lb) where the 
simulation starts with a pure mode with £ = m = 1 and Ur = 0. The initial amplitude 
of the hrst peak is chosen to be 0.40 and 0.45 in the cases (la) and (lb), respectively. 

In the simulation (la), the nonlinear effect is not so strong. Figure shows the 
prohle of the scalar held on the 0 = 0 line in the equatorial plane at hve moments when 
the scalar held peak passes the 0 = 0 line. During 0 < f/M < 1000, the axion cloud 
becomes concentrated around the black hole and the peak value becomes larger. But 
after that, the axion cloud recedes from the black hole, and the prohle at t/M ^ 2000M 
resembles the initial prohle. Therefore, in this case, the nonlinear self-interaction ehect 
just causes the radial oscillation of the axion cloud with the period of ~ 2000M. 
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t= \1M, 507M, 999M, 1490M, 198IM 



rJM 


Figure 6. Snapshots of the scalar field on the ^ = 0 line in the equatorial plane for 
t = 17M (solid line, red online), 507M (long-dashed line, green online), 999M (short- 
dashed line, blue online), 1490M (dotted line, purple online), and 1981M (dash-dotted 
line, sky blue online) in the simulation (la). Time is chosen in order to depict the 
moments when the scalar field peak passes the (j) = 0 line. 


The superradiant instability continues throughout the simulation (la), and the 
averaged rate of the energy extraction is given by {dE/dt) 4 x 10“^(/a/Mpi)^ with 
the Planck mass Mpi ~ 10^® GeV. This value is larger by a factor of ~ 3 than that 
of the Klein-Gordon held with the same amplitude due to the nonlinear self-interaction 
effect. The energy extraction is expected to change the axion cloud energy from that 
of (la) to that of (lb), which are E/M ^ 2.6 x 10^(/a/Mpi)^ and 3.3 x 10^(/a/Mpi)^, 
respectively. The required time for the growth from (la) to (lb) is estimated to be 
< 1.7 X lO^M. 

In the simulation (lb), the held behaviour is fairly diherent from the case of (la). 
Figure [^shows the snapshots of the scalar held value in the equatorial plane at t = 0 (top 
left), 800M (top right), lOOOM (middle left), 1150M (middle right), 1400M (bottom 
left), and 1900M (bottom right). The axion cloud gradually becomes concentrated 
around the black hole {t = 800M), and around t = lOOOM, the axion cloud becomes 
strongly distorted by the excitation of various modes, especially the m = — 1 mode that 
falls into the black hole. After that, the axion cloud recedes from the black hole for a 
while {t = 1150M). The axion cloud in this intermediate state consists of two compact 
solitonlike objects slowly rotating around the black hole. These solitons approach the 
black hole again around t = 1400M, exciting various modes. After that the axion cloud 
settles down to a state somewhat similar to the initial state but still continuing the 
excitation of the m = — 1 mode {t = 1900M). This is the bosenova collapse in the case 
Mfi = 0.30. 

Figure shows the energy and angular momentum fluxes towards the horizon 
observed at r* = —200M as functions of time. The inset enlarges the period 0 < t/M < 





Bosenova and Axiverse 


26 


t/M = 000.0 


t/M = 800.0 
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Figure 7. Snapshots of the scalar field in the strongly nonlinear, £ = m = 1 case 
[simulation (lb)] at / = 0 (top left), 800M (top right), lOOOM (middle left), 1150M 
(middle right), 1400M (bottom left), and 1900M (bottom right) in the equatorial plane 
in the simulation (lb). 


1000. In this period, both the energy and angular momentum fluxes are negative, and 
the energy is extracted from the black hole. During 1000 < t/M < 2000, there are 
two periods when large positive energy rapidly falls into the black hole, corresponding 
to the high concentration of the axion cloud at t/M k, 1000 and 1400. The bosenova 
collapse is characterized by the excitation of the m = — 1 mode, which carries positive 
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Figure 8. Energy flux (solid line, red online) and angular momentum flux (dashed 
line, green online) towards the horizon observed at r* = —200M as functions of time 
in the simulation (lb). The inset enlarges the period 0 < t/M < 1000. 


energy and negative angular momentum to the black hole terminating the superradiant 
instability. During this period, about 5.9% of the total energy falls into the black hole. 


/ = 1950M 



rJM 


Figure 9. Energy density (solid line, red online) and angular momentum density 
(dashed line, green online) with respect to the tortoise coordinate r* in the distant 
region at t = 1950M in the simulation (lb). 

Next, let us look at the scattering of the scalar held to the distant region. Before the 
bosenova, the held conhguration is practically conhned in the region r* < 200M, and 
it is approximately zero beyond this region due to the exponential decay. But after the 
bosenova, the axion cloud begins to spread out to the distant domain. Figure shows 
the energy and angular momentum densities with respect to the tortoise coordinate r* 
in the distant region at t = 1950M. About 13.4% of the total energy distributes in the 
domain r* > 300M. 
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r= 1950M 


t= 1950M 




rJM 


rJM 


Figure 10. Left panel: Snapshots of the real part of the scalar field components 
for m = 1 and n = 0 (solid line, red online), 2 (long-dashed line, green online), 4 
(short-dashed line, blue online), and 6 (dotted line, purple online) as functions of the 
tortoise coordinate r* in the distant region at t = 1950M in the simulation (lb). Right 
panel: The same as the left panel but for to = 1 (solid line, red online), 3 (long-dashed 
line, green online), 5 (short-dashed line, blue online), and 7 (dotted line, purple online) 
and n = 0. 


Whether these scattered fields to the far region are gravitationally bounded or not 
is of interest because it affects the evolution of the axion cloud in a longer time scale. 
The exact separation of bounded and unbounded modes requires mode decomposition 
both in time and radial directions, and we have not calculated it. But we can understand 
the gross feature by looking at the waveform because bounded modes have frequencies 
smaller than the mass /i and decay exponentially at large r*, while unbounded modes 
oscillate with frequency uj larger than /i and decay slowly. Figure IT shows snapshots 
of for m = 1 and n = 0, 2,4, 6 (left panel) and for m = 1, 3, 5, 7 and n = 0 (right 
panel) at t = 1950M as functions of the tortoise coordinate r*. Clearly the major 
component of the m = 1 mode oscillates with frequencies larger than /i. It is also found 
that the modes with m > 3 include frequencies larger than /r, although some waves have 
frequencies close to /i. Therefore, fairly large part of the fields scattered to the distant 
region is unbounded, and they propagate towards infinity. This is in contrast to the 
case of Mfi = 0.40 studied in na, where only a small portion of the field is scattered 
to infinity and most part of the scattered fields remains bounded. In this sense, the 
bosenova explosion becomes more violent as Mp is decreased. 

In a more realistic situation, the bosenova is expected to happen when the energy 
amount close the case (la). From this, the condition for the bosenova occurrence is 
estimated to be Ea/M ^ 3 x 10^(/a/Mpi)^ for Mp = 0.30. This energy amount is 
larger than that for Mjj, = 0.40 by a factor of two. The axion cloud energy at the 
bosenova occurrence becomes larger as M/i is decreased. This is consistent with the 
naive expectation that the bosenova happens when the peak field amplitude becomes of 
the order /a, because in that case, we obtain the estimate E/M ~ fx^ip'^V/M oc 
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4-2.2. Pure i = m = 2 mode For a* = 0.99, the fundamental i = m = 2 mode has 
the fastest growth rate for 0.448 < Mfj, < 0.928. In this section, we choose Mfj, = 0.80 
and simulate the evolution of the axion cloud starting from the fundamental £ = m = 2 
mode whose first peak value is 1.00. 
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Figure 11. Snapshots of the scalar field value in the equatorial plane in the £ = m = 2 
case [simulation (2)] at t = 0 (left) and lOOOM (right). 

Figure [TT| shows the snapshots of the field value on the equatorial plane at t/M = 0 
and 1000. Because we consider the £ = m = 2 fundamental mode, there are two 
local maxima and two local minima. Due to the nonlinear interaction, the axion cloud 
becomes closer to the black hole during the evolution. This is a kind of relaxation 
process, and the axion cloud remains concentrated around the black hole after that. 
The local maxima keep fairly large values that are between two and three. The typical 


conhguration is depicted in the right panel of figure This phenomenon is in contrast 
to the £ = m = 1 case. In the simulation (la) in the previous section, the axion cloud 
starting from the pure £ = m = 1 mode oscillates in the radial direction. After a 
mild concentration near the black hole, the axion cloud recedes from the black hole, 
and the configuration at f ~ 2000M is very close to the initial configuration. In the 
£ = m = 2 case, the axion cloud never resembles the initial configuration and keeps high 
concentration around the black hole. 

Figure [T^ shows the energy and angular momentum fluxes towards the horizon 
observed at r*/M = —100. The fluxes during 100 < t/M < 300 are generated by 
the relaxation process. After this relaxation process, the energy flux is always negative 
and the energy extraction continues. Periodically, the absolute value of the energy flux 
becomes large. This is because the axion cloud shows a small radial oscillation. When 
the axion cloud becomes closer to the black hole, the amplitude of the field oscillation 
increases to ~ 3, and at that moment large amount of energy flux falls into the black 
hole. When the axion cloud becomes a bit far from the black hole, the amplitude of the 
field oscillation decreases to ~ 2 and the absolute value of the energy flux also becomes 
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Figure 12. Energy flux (solid line, red online) and angular momentum flux (dashed 
line, green online) towards the horizon observed at r* = — lOOM as functions of time 
in the simulation (2). The inset enlarges the period 500 < t/M < 2000. 

small. 


t = 2000M 



rJM 


Figure 13. Energy density (solid line, red online) and angular momentum density 
(dashed line, green online) with respect to the tortoise coordinate r* in the far region 
50 < r*/M < 2000 at t = 2000M in the simulation (2). 


This behaviour of the energy flux towards the horizon is analogous to the weakly 
nonlinear case in the ^ = m = 1 case. But an important difference appears in the held 
behaviour at the distant region. Figure 13 shows a snapshot of the energy and angular 
momentum densities at t/M = 2000 at far region, 50 < r*/M < 2000. About 15% 
of the total energy distributes in this domain. This is in contrast to the i = m = 1 
case, where the amount of the energy scattered to the distant region is very small (see 
hgure 7 of na). Therefore, in spite of the fact that the bosenova collapse (i.e., the 
positive energy infall towards the horizon) does not happen, the scattering of the held 
to the distant region occurs continuously. As a result, the axion cloud in the domain 
r^/M < 50 continuously loses its energy. 
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t/M 

Figure 14. The increse in the energy in the distant region Ai^dist (solid line, red 
online) and the extracted energy (dashed line, green online) from the black hole AEg^t 
from time lOOOM to t in the simulation (2). See text for details. 


To see this fact more clearly, let us define the following quantities: 


A-E'dist -£^dist(^) ~ -Sdist(lOOOM), 
rt 




\Fe\ dt. 


(4.9a) 

(4.96) 


' lOOOM 


Here, T^dist(^) denotes the energy amount in the domain r* > 50M, and therefore, AT^dist 
indicates the energy amount that is supplied from the central region r* < SOM between 
time lOOOM and t. On the other hand, Ai^ext indicates the amount of the extracted 


energy from the black hole between lOOOM and t. Figure shows the behaviour of 
these two quantities as functions of time. We see that although some fluctuation exists, 
the distant energy has tendency to increase in time, and its amount is much larger than 
the extracted energy. 

The important question is whether the scattered held at the distant region is 
gravitationally bounded or not. In order to derive the answer, one must decompose 
the modes in both time and radial directions, and we have not calculated it. Here, 


we discuss the gross feature from hgures 14 and 13 In hgure the value of AF^dist 


sometimes decreases in time, and this means that the scattered helds include bounded 
modes. But since the gross tendency is that AF^dist increases in time, the scattered helds 


also include the nonnegligible fraction of unbounded modes. Figure also indicates 
that some energy generated after the relaxation process propagates towards the very far 
region. 

Let us discuss the implication of this simulation. Since the main part of the axion 
cloud loses its energy while the energy extraction from the black hole continues in this 
simulation, the setup here is unrealistic in the sense that it does not approximate the 
end point of the superradiant instability. In a realistic situation, when the amount 
of the axion cloud energy is small, the nonlinear ehect is not important and the 
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axion cloud grows by the superradiant instability. As the nonlinear effect becomes 
important, the scattering towards the far region gradually occurs, while the extraction 
of the energy from the black hole continues. At some energy, the scattering rate of 
energy to the distant region and the rate of the energy extraction from the black hole 
would balance, and the final state would not be the bosenova collapse, but rather a 
steady outflow of the field from the axion cloud. This steady state should be realized 
before the axion cloud energy reaches the energy amount in the simulation here, that 
is, E/M ^ 1.9 X 103(/,/Mpi)2. 

Note that we have performed several simulations by changing the initial amplitude. 
If we choose an artificially large initial amplitude such as 1.2, the phenomena analogous 
to the bosenova in the i = m = 1 case occurs. There, the excitation of the m = —2 
mode is observed and it continues for the period of ~ 400M. However, we consider 
this “bosenova collapse” would not happen in a realistic case, because the axion cloud 
energy cannot reach the required energy to cause the bosenova due to the energy loss 
by the scattering of the field to the distant region. 


4-2.3. Superposition of i = m = 1 and 2 modes In this section, we study the case 
that the axion cloud consists of the superposition of the fundamental i = m = 1 and 
i = m = 2 modes initially. Here, the primary component of the axion cloud is the 
i = m = 1 mode, and the i = m = 2 mode is added as a perturbation. To be more 
specific, we choose the case of a* = 0.99 and M/i = 0.40, and assume that the initial 
peak value of the i = m = 1 mode to be 0.70. If the £ = m = 2 mode is absent, the 
system is same as the one simulated in our previous paper [T^. To this system, we 
add the £ = m = 2 mode whose first peak value is 0.01. Our primary interest here is 
what happens to the £ = m = 2 mode under the influence of the nonlinear effect of the 
£ = m = 1 mode. To see this, it is useful to calculate the decomposition 


^ ^ imjt ^ r A)YiYn{9 ^Im" 


(4.10) 


im 


Then, 622 approximately gives the £ = m = 2 component of the axion cloud. 

Figure lA shows snapshots of the real part of 2622 F 22 on the 0 = 0 line in the 
equatorial plane at t/M = 0, 399, and 605. Although the added £ = m = 2 mode is 
small initially, during the bosenova, it grows to become order one analogously to the 
resonance by a forced oscillation {t = 399M). The profile of the £ = m = 2 mode 
becomes highly concentrated around the black hole, and generated ingoing waves carry 
positive energy flux to the horizon. As the bosenova ceases, the £ = m = 2 mode 
also settles down, but the mode function is very different compared to the initial state 
{t = 605M). The £ = m = 2 mode after the bosenova has a larger frequency compared 
to the initial state and spreads over a distant region, causing several oscillations in the 
radial profile of the mode function. In this phase, infalling waves to the horizon include 
both the m = ±2 modes, and the net flux carried by the £ = m = 2 mode remains 
positive at least up to t = lOOOM. 
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f/M=0, 399, 605 



Figure 15. Snapshots of the real part of the component 2bimYim for £ = m = 2 on 
the (j> = 0 line in the equatorial plane at t/M = 0 (solid line, red online), 399 (long- 
dashed line, green online), and 605 (short-dashed line, blue online) in the simulation 
(3). At t/M = 399, the field is concentrated around the black hole and very high peak 
value. After that, the field settles to a configuration that resembles an overtone mode 
{t/M = 605). 


As observed above, the evolution of the £ = m = 2 mode is quite complicated. 
Although the superradiant instability by the i = m = 2 fundamental mode no longer 
occurs after the bosenova, the m = 2 mode survives showing complicated dynamics. 
With the simulation results up to t = lOOOM, it is very difficult to predict the hnal fate 
of the i = m = 2 mode in a long-term evolution. 


t/M =503 



rJM 

Figure 16. Snapshots of the scalar field $ on the b = 0 line in the equatorial plane 
at t = 503M in the cases that the small ^ = m = 2 is present (solid line, red online) 
[simulation (3)] and absent (dashed line, green online) [15]. 
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Another interesting phenomena observed in the simulation is that the property of 
the mode excitation during the bosenova is rather changed by the presence of the small 
i = m = 2 mode. Figure 1^ shows a snapshot of the held value <F on the line 0 = 0 
in the equatorial plane at t = 503M. The cases that the i = m = 2 mode is present 
and absent are shown by the solid line and the dashed line, respectively. When the 
small i = m = 2 mode is added, the excited mode has larger amplitudes and the typical 
frequency becomes smaller. 




t/M HM 

Figure 17. The energy flux Fe (left panel) and the angular momentum flux Fj 
(right panel) towards the horizon observed at r* = —200M as functions of time for 
the cases that the small i = m = 2 \s present [simulation (3)] (solid line, indicated by 
“Superposed case”) and absent [TS] (dashed line, indicated by “Pure case”). 

Corresponding to this phenomena, the energy and angular momentum fluxes 
towards the horizon also change as shown in hgure In this hgure, we compare the 
energy flux (left panel) and the angular momentum flux (right panel) towards the horizon 
observed at r*/M = —200 in the two cases that a small £ = m = 2 mode is present (solid 
line) and absent (dashed line). The energy flux becomes much larger while the absolute 
value of the angular momentum flux becomes smaller when the small £ = m = 2 mode 
is added. Therefore, the evolution of the energy and angular momentum of the axion 
cloud is strongly affected by the mixture of additional modes to the axion cloud even if 
they are small. 


4-2.4- Superposition of = 0 and 1 modes Here, we consider the case in which the 
initial state is given by the superposition of the fundamental mode (ur = 0) and the 
overtone mode {n,. = 1) for £ = m = 1. As the system parameters, we choose a* = 0.99 
and Mfi = 0.40, and prepare the initial data by adding the two modes = 0 and n,- = 1 
with the hrst peak values 0.60 and 0.20, respectively. When only the fundamental mode 
rir = 0 is present with the peak value 0.60, the bosenova does not happen [TB]. We are 
interested in what is changed if the n,- = 1 mode is superposed to this situation. 

Figure 18 shows the held conhguration in the equatorial plane at t/M = 0, 400, 
and 700. In our initial data, the inner local maximum (minimum) exists at r^/M ^ 12.5 
and 0 = 0 (tt) and the outer local minimum (maximum) exists at r^/M k. 77.5 and 
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Figure 18. Snapshots of the scalar held in the case that the two modes nj- = 0 and 
rir = 1 are superposed [simulation (4)] at t = 0 (top left), 400M (top right) and 700M 
(bottom) in the equatorial plane. 


0 = 0 (vr) in the equatorial plane. During the evolution, it was found that the outer 
local maximum/minimum outside is absorbed to the inner local maximum/minimum 
{t/M = 400). In other words, the merger of the two maxima/minima happens. Then, 
the system resembles an axion cloud with just one fundamental mode {t/M = 700). 
Therefore, a kind of conversion from the overtone mode to the fundamental mode occurs 
due to the nonlinear self-interaction. 

The evolution after this merger is very analogous to the case of pure fundamental 
mode with i = m = 1 with initial peak value 0.70 simulated in |T5]. After high 
concentration around the black hole, various modes are excited. In particular, the 
m = — 1 mode is the primary mode that falls into the horizon. Figure depicts 
the held conhguration at t = 1210M on the 0 = 0 line in the equatorial plane. The 
excitation of the m = — 1 mode is clearly observed. 

To summarize, the existence of the overtone mode with Ur = 1 does not much 
affect the dynamics of the axion cloud which is dominated by the fundamental mode 
with Ur = 0, because the overtone mode is converted to the fundamental mode due to 
the nonlinear self-interaction effect. Our results suggest that when the nonlinear effect 
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t=moM 
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Figure 19. A snapshot of the scalar field in the case that the two modes rir = 0 and 
rir = 1 are superposed [simulation (4)] as a function of r^/M at t = 1210M on the 
(j) = 0 line in the equatorial plane. 


is strong, it would be sufficient to consider the dynamics of the fundamental mode. 

4-3. Summary of this section 

In this section, we numerically studied the phenomena caused by the nonlinear 
interaction for various setups. For the i = m = 1 mode, we studied the dependence 
of the bosenova phenomena on the value of Mfj,. In the case of Mjx = 0.30 simulated 
in (lb), the axion bosenova happens when the axion cloud acquires a larger energy 
compared to the case of Mfi = 0.40 studied in [16]. The bosenova phenomena becomes 
more violent for a smaller M/i, causing stronger excitation of gravitationally unbounded 
modes. The excitation of the m = —1 mode is commonly observed for both M/i = 0.30 
and 0.40. The m = —1 mode falls into the black hole carrying positive energy, and thus, 
it terminates the superradiant instability. 

In the i = m = 2 case, we found the evidence against the bosenova occurrence. 
In the simulation (2), no phenomena analogous to the bosenova happens, and instead, 
we observed the continuous energy extraction from the black hole and generation of a 
steady outflow to the distant region, simultaneously. For a chosen parameter in the 
simulation (2), the outflow of the scalar held energy was larger than the rate of the 
energy extraction. From these results, we suspect that the growth of a realistic axion 
cloud formed by the superradiant instability would stop when the energy extraction 
rate and energy outhow rate balance with each other. Note that if we choose a large 
initial amplitude by hand, the phenomena analogous to the bosenova was observed. In 
particular, we observed the excitation of the m = —2 mode, which is an analogous to 
the m = — 1 mode excitation in the i = m = 1 case. However, in a realistic situation, 
the axion cloud would not be able to acquire the required energy to cause a bosenova 
collapse due to the generation of a steady outflow. 
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In the simulations (3) and (4), we studied the cases where two modes are 
superposed. In the simulation (3), we added the £ = m = 2 fundamental mode as 
a perturbation to the axion cloud in the i = m = 1 mode. During the bosenova, the 
i = m = 2 mode also carries a positive energy flux to the horizon, but the forced 
oscillation by the bosenova increases the energy of the i = m = 2 mode. The mode 
function after the bosenova is very different from the initial state, and spreads over a 
larger region. It was also observed that the added i = m = 2 mode strongly affects 
the properties of the bosenova, such as the values of the energy/angular momentum 
fluxes towards the horizon. On the other hand, the simulation (4) indicates that the 
addition of the overtone mode = 1 to the fundamental = 0 mode does not strongly 
affect the properties of the bosenova, because the overtone mode is converted into the 
fundamental mode by the nonlinear self-interaction. 

5. Gravitational wave emission 

One of the important phenomena caused by an axion cloud around a black hole 
is the gravitational wave emission. In our previous research, we calculated the 
amplitudes of continuous gravitational waves emitted by a Klein-Gordon scalar cloud 
ignoring the nonlinear self-interaction [16], and discussed observational constraints on 
axion model parameters im. The neglect of the nonlinear self-interaction in this 
calculation may, however, result in an estimation largely different from reality, because 
the axion bosenova may emit burst-type gravitational waves with amplitudes much 
larger than those of continuous emissions, and besides, the self-interaction may cause 
fluctuation/modulation in the amplitude/frequency of continuous waves emitted in the 
superradiant phase, which has significant influences on the wave search in the data 
analysis. Thus, more elaborate calculations of gravitational waves taking account of the 
nonlinear self-interaction are highly desired. In this section, we describe the current 
status of our research in this direction. 


5.1. Formulation and numerical method 


We first explain the basic strategy of the calculation. In the previous section, we 
have numerically determined the behaviour of the axion scalar field in the test field 
approximation. Using this numerical solution, the energy-momentum tensor can be 


calculated using (3.12). Then, the amplitudes of gravitational waves can be determined 


by solving the linearized Einstein equations with the source term in the Kerr 
background. Throughout the present paper, we ignore the radiation reaction to the 
dynamics of the scalar field and the black hole parameters. 

Various approaches can be taken to determine the metric perturbation for a given 
source term. In the present paper, we choose the approach to directly solve the Teukolsky 
equation [40] in the time domain. In this approach, the perturbed Einstein equations 
in the Kerr background are reduced to a set of PDEs for the perturbation of the 
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spin coefficients and the Newman-Penrose quantities i/'o and "04 constructed from the 
Weyl curvature tensor. This set of PDEs can be further reduced to decoupled second- 
order master wave equations for two gauge-invariant complex functions ± 2 ^^ that are 
algebraically related to 0o and 04 , known as the Teukolsky equations, and the relation 
between the two Teukolsky functions ± 2 ^- With the help of this relation and the original 
PDEs, we can completely determine all Newman-Penrose quantities as well as the metric 
perturbation from either of these two functions 

The explicit form of the Teukolsky equation for a gravitational perturbation m is 
given by 
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+ (s2cot2 0-s),T = 47rST (5.1) 


with s = +2 or —2. The right-hand side is the source term, and T is given in terms of 
components of the energy-momentum tensor that generates the gravitational waves. In 
our case, we choose s = —2 for a technical reason in numerical calculations. For this 
choice of s. 


T = (2Ip‘)T,, 


(5.2) 


with 


T4 — (^ + 37 — 7* -|- 4/i -|- — 2t* + 2a)Tnm* — (A -|- 27 — 27* -|- 

+(^* — r* -f- 0* -|- 3 q; -|- 47r)[(A -|- 27 -|- 2fi*)Tnm* ~ — t* + 2/3* + 2q;)T„„], (5-3) 


where T^n 


T^ipU^n'^, and so on. Here, and are the Kinnersley null tetrads. 
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(5.4a) 

(5.40 

(5.4c) 


A and 5 are the derivative operators, A := and 5 := m^V^, and a, 0 , 7 , /i, and 

TT are non-vanishing Newman-Penrose variables for the spin coefficients m 

The Teukolsky equation has terms with divergent coefficients proportional to 
l/sin^0, which may cause a problem in numerical simulations. In order to resolve 
this, we spectrally decompose _ 24 / as 

_2T= ^ —0('^)(f,r,0)sinl™-2|-cosl™+2|(5.5) 

m=o(^±2,.. ^ 2 2 


Note that the this expression is consistent with the behaviour of the spin-weighted 
spheroidal harmonics if 0 ('") is hnite and regular at the two poles. The extra factor 
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A^/r is chosen so that the amplitudes of for ingoing waves to the horizon and for 
outgoing waves to inhnity become constant. If the source term is similarly decomposed 
as 

T = (5.6) 

m=0,±l,±2,... 

we have the equation 
- [(r^ + ay - Aa\l - + (r^ + 

rA 


—4 [My — a?) — A(r + iay) + imMar^^ + 2 (r^ + a^) 
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The left-hand side of this equation is explicitly regular. Instead, the divergent coefficient 
appears in the right-hand side. The regularity of this source term is discussed below. 

Now, let us discuss the source term in more detail. It is convenient to adopt the 
decay constant /a as the unit of the scalar held $ as done in (4.1), and we use = <h//a 
below. This means that the Teukolsky variable is normalized by f^. For the energy- 


momentum tensor of the scalar held (3.12), each component in ( |5.3 ) becomes 
T„m- = r„.„. = T„„ = 


(5.8) 


The important point is that since and are null and orthogonal to each other, 
there is no contribution from the potential term U{(p). We substitute (4.4) into (5.8) 
and then, into (5.3). Since T 4 is quadratic in the scalar held (p, we label one of this by 
mi and the other by m 2 . Then, T can be expressed in the form 


T = 


mi ,m2)p,i{mi+m2)(t> 


(5.9) 


mi,m2 


This relation implies that the coupling between the mi mode and the m 2 mode of the 
scalar held generates gravitational waves in the m = mi-|-m 2 mode, denoted by 
Comparing this expression with (5.6), we hnd that gravitational waves in the m mode 
can be expressed as the sum of the contributions from all pairs of modes (mi, m 2 ) 
satisfying the condition m = mi -|- m 2 . 

We have to note the subtlety in the meaning of the mode number m. As noted in 
the previous section, the mode function of the scalar held in (4.4) includes both 
positive and negative frequency modes, which can be reinterpreted as the true -|-m and 
—m modes, respectively. As a result, the Teukolsky function contains 

both positive and negative frequency modes, which are naturally reinterpreted as the 
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+rh and —fh modes, respectively. Therefore, in order to specify the true +m mode, we 
have to decompose positive and negative frequency modes in 

and superpose the positive frequency mode in negative frequency mode 

in In this paper, we do not calculate this decomposition and resummation, 

because the negative m mode of the scalar field gives a minor contribution to the axion 
cloud in the Schwarzschild case considered below. But this procedure would be required 
in the more realistic Kerr case. 


Since the denominator of the right-hand side of (5.7) becomes zero in general at 


y = ±1, its regularity has to be checked and a regular expression must be implemented 
in the numerical code. In the Kerr case, we checked the regularity at the poles using 
MATHEMATICA, although we do not present the explicit results here because it is very 
tedious. In the Schwarzschild case, we can find a relatively simple expression, and it is 
presented in Appendix B[ 

Similarly to the pseudospectral approach in the scalar case, we further decompose 
as 




n=0,l,2,... 


a, 


(m) 




( 6 . 10 ) 


By substituting this expression into the left-hand side of (5.7), and the corresponding 
series expressions for and respect y, (4.6), into the right-hand side. 


and by equating the coefficients of each y"', we obtain a sequence of (1 -|- l)-dimensional 
PDEs for with sources expressed in terms of known functions of More detailed 


formulas are shown in Appendix B[ 

In simulations, we first calculate the functions for the scalar field. Then, for 
fixed mi and m 2 , we solve the functions by using the data and 

for the source term. We use the sixth-order hnite differencing in spatial direction and 
the fourth-order Runge-Kutta method in time direction. In this method, a numerical 
instability develops around 0 < r* < lOOM. In order to stabilize the numerical data. 


we added a stabilization term to the left-hand side of (5.7), where 


5(r*) = ei exp 




Wi 


€2 exp 


42 )' 


W2 


(5.11a) 


with 


ei = 0.015, 62 = 0.005, = 20M, rf) = 70M, wi = W 2 = SOM. 


(5.116) 


Because adding such a term is artihcial, there remains a room for improvement. Here, 
we postpone this issue as a future research, and show the results obtained under this 
method. But we would like to point out that because generated gravitational waves 
propagates towards the horizon or infinity, the period during which the stabilization 
term affects the gravitational wave form is relatively small. Note that in the scalar case 
in the previous section, it was not necessary to add such a stabilization term. This 
problem is specific to the Teukolsky equation for a gravitational perturbation. 
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5.2. Numerical results for the Schwarzschild case 

At this point, we have not yet developed a code for a Kerr black hole case because 
the source term is very complicated. In this paper, as a preliminary research, we 
present the results for the Schwarzschild case. Because the energy extraction does not 
happen, we cannot discuss the generation of gravitational waves at bosenova as a result 
of the superradiant instability. However, in the case of Klein-Gordon field without self¬ 
interaction in the Schwarzschild background, there exist long-lived quasibound states 
with small decaying rates: For M/i = 0.30, the decay rate of the i. = m = 1 fundamental 
mode is Mui ~ —0.9456 x 10“^. Then, we can prepare a Klein-Gordon quasibound state 
with a certain amplitude as the initial condition by hand, and evolve the scalar field by 
switching on the nonlinear self-interaction. The phenomena in this setup captures some 
of the features of the axion bosenova in the more realistic Kerr-background case within 
the time scale of few thousand M, and therefore, we can understand the gross feature 
of gravitational waves emitted during the bosenova. 

In what follows, we show the numerical results for the three cases: the Klein-Gordon 
case (without nonlinear self-interaction), the mildly nonlinear case, and the strongly 
nonlinear case, one by one. The mass of the scalar field is fixed to be M/i = 0.30, and 
we use the scalar field in the quasibound state in the i = m = 1 fundamental mode as 
the source for gravitational waves. 

Note that if the black hole is a nonrotating Schwarzschild black hole and the scalar 
field is symmetric with respect to the equatorial plane as assumed above, is 

related to as r,//) = r, —//). Therefore, it is sufficient to consider 

the positive number of rh. 


5.2.1. Klein-Gordon case We first briefiy discuss gravitational wave emissions from the 
Klein-Gordon quasibound state. In this case, the spectral component for gravitational 
waves, with {n,m) = (0,2), is generated from the scalar field with {n,m) = 
(0, ±1). No excitation of higher n modes has been found, which is consistent with our 
previous perturbative study [T6] . 

shows snapshots of the real and imaginary parts of with {n,rh) = 


Figure 


20 


(0, 2) for t = lOOM (Left panel) and t = 400M (Right panel). We started with the initial 
data = 0. Since they are not realistic initial data in the sense that they do 

not reflect continuous waves emitted before f = 0, a spurious gravitational wave burst 
can be seen in the initial stage. But they are soon radiated away, and the generation of 
continuous waves are observed after a sufficient time. 


5.2.2. Mildly nonlinear case Now, we consider the case in which the scalar field follows 
the sine-Gordon equation. Here, we consider the mildly nonlinear case, where the initial 
amplitude of the scalar field oscillation is chosen as (ppeak = 0.40. 

Before looking at generated gravitational waves, let us inspect the scalar field 
behaviour. Figure shows the snapshots of the scalar field on the line 0 = 0 in 
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t = lOOM t = 400M 




Figure 20. Snapshots of the real part (solid line, red online) and the imaginary part 
(long dashed line, green online) for (^("“^ for m = 2 in the equatorial plane) at 
time t = lOOM (Left panel) and t = 400M (Right panel) in the Klein-Gordon case. 


t = AM, 493M, 1003M, 1496M 



rJM 


Figure 21. A snapshot of the scalar field on the (j> = 0 line in the equatorial plane for 
t = AM (solid line, red online), A93M (long-dashed line, green online), 1003M (short- 
dashed line, blue online), and 1496M (dotted line, purple online). Time is chosen in 
order to depict the moment when the scalar field peak passes the (p = 0 line. 


the equatorial plane at time t = AM (solid line), 493M (long-dashed line), 1003M 
(short-dashed line), and 1496M (dotted line). Here, we selected the moments at which 
the scalar held peak passes the 0 = 0 line. Although the held amplitude changes by 
a factor due to self-interaction, the system does not experience a violent phenomena 
analogous to the bosenova. In fact, the amplitude increases up to t ~ lOOOM, but it 
decreases after that. As a result, the scalar held distribution is very similar for t = A93M 
and 1496M. Therefore, the situation is similar to the weakly nonlinear case discussed 
in the rotating background in [13]. In the pseudospectral approach, the held primarily 
consists of with (n, m) = (0,1), and other components are scarcely excited. 

Although the behaviour of the scalar held caused by the nonlinear self-interaction is 
not very dynamic, interesting phenomena appear in the generated gravitational waves. 
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Below, we show just the gravitational wave forms generated by mi = m 2 = ±1 mode of 
the scalar held, because the higher m modes remain to be small. Therefore, we consider 
the gravitational waves in the m = 2 mode excited by the mi = m 2 = 1 mode of the 
scalar held. In the mildly nonlinear case, gravitational waves in the m = 0 mode excited 
by the mi = —m 2 = ±1 modes are scarcely observed. 
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Figure 22. Snapshots of the real (solid line, red online) and imaginary (long dashed 
line, green online) part of spectral component Oq ' of the Teukolsky variable for a 
gravitational perturbation at time t = 600M (top left panel), t = 950M (top right 
panel), t = 1150M (bottom left panel) and t = 1500M (bottom right panel) in the 
mildly nonlinear case. 


Figure 22 shows snapshots of the real and imaginary parts of with (n,m) = 


(0,2) for t = 600M (top left panel), t = 950M (top right panel), t = 1150M (bottom 
left panel), and t = 1500M (bottom right panel). Note that this quantity agrees with 
the value of in the equatorial plane (t/ = 0). Similarly to the Klein-Gordon case, 
continuous waves appear after the radiation of spurious initial burst. But in this case, 
its amplitude increases in time. Around t = 900M, the “beating” pattern appears in the 
radiated gravitational wave, and it becomes more prominent as time goes. The “beat 
frequency” changes in time and becomes smaller as time goes, as can be understood by 
comparing thn three panels for t = 950M, 1150M, and 1500M. 

(h, m) = (0, 2) as a function of time observed 


Figure 


23 


shows the value of for 


at the hxed position r* = 200M. After a linear growth in the amplitude, the beating 
phenomena are observed. 

Figure 


24 


shows a snapshot for t = 1500M, but for with m = 2 and n = 0 
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= 200M 



( 2 ) 

Figure 23. Behaviour of real part of the spectral component ap of the Teukolsky 
variable for a gravitational perturbation as a function of time observed at the point 
r* = 200M in the mildly nonlinear case. 


t= 1500M 



rJM 


Figure 24. A snapshot of the real part of the spectral component of Teukolsky 
variable for a gravitational perturbation for fixed m = 2 and n = 0 (solid line, red 
online), 1 (long-dashed line, green online), 2 (short-dashed line, blue online), and 3 
(dotted line, purple online) at time t = 1500M in the mildly nonlinear case. 


(solid line), 1 (long-dashed line), 2 (short-dashed line), and 3 (dotted line). In contrast 
to the case of the Klein-Gordon field, the higher n modes are also excited. This means 
that gravitational waves in the modes with i > 2 are excited due to nonlinear self¬ 
interaction of the scalar held when they are expanded by the spin-weighted spheroidal 
harmonics 

To summarize, in the mildly nonlinear case, several modihcations from the Klein- 
Gordon case appear in the gravitational wave forms although the scalar held is not 
highly dynamical. After a growth in the amplitude of monochromatic waves, beating 
pattern appears at this point. The beating frequency changes in time: In our simulation 
time scale, it becomes smaller as t is increased. Also, gravitational waves in the higher 
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i modes are excited. 
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Figure 25. Snapshots of the scalar field in the strongly nonlinear case at / = 0 
(top left), 600M (top right), 750M (bottom left), and 850M (bottom right) in the 
equatorial plane. 


5.2.3. Strongly nonlinear case Now, we turn our attention to the case of the strongly 
nonlinear case. Here, we choose the initial amplitude of the scalar held oscillation to be 
‘^peak = 0.5. Figure 26 shows the scalar held conhguration in the equatorial plane. The 
scalar held becomes highly concentrated towards the black hole due to nonlinear self¬ 
interaction rapidly, and after several turns, the concentration is destroyed, and the held 
conhguration changes to a disperse spiral pattern. This phenomenon is quite analogous 
to the bosenova observed in the rotating case in our previous work: Compare this hgure 
with hgure 8 of na. Figure depicts twice of the real part of the spectral component 
of the scalar held with m = 1 and n = 0 (solid line), 2 (long-dashed line), 4 
(short-dashed line), and 6 (dotted line). In this case, the scalar held in the m = 1 
gets amplihed and falls into the black hole, and no excitation of the m = — 1 mode is 
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t = 825M 



rJM 


Figure 26. A snapshot of the twice of the real part of the spectral components 
for the scalar field with m = 1 and n = 0 (solid line, red online), 2 (long-dashed line, 
green online), 4 (short-dashed line, blue online), and 6 (dotted line, purple online) at 
time t = 825M. This depicts the moment at which a large amount of the scalar field 
energy falls into the black hole. 


observed. This is a slight difference from the realistic bosenova in the rotating black 
hole case [T3] . 

Now, we study gravitational waves. Similarly to the mildly nonlinear case, we study 
just the gravitational wave forms generated by the rui = m 2 = ±1 modes of the scalar 
held. But in this case, we consider both the rh = 2 mode excited by the mi = m 2 = 1 
mode and the m = 0 mode excited by the mi = —m 2 = ±1 modes. Figure ^7\ shows 


snapshots of the real and imaginary parts of with (h, m) = (0, 2) for t = 400M (top 
left panel), t = 650M (top right panel), t = 850M (bottom left panel), and t = 1200M 
(bottom right panel). Because the gravitational wave amplitude changes signihcantly in 
time, we show these quantities in the logarithmic scale. Similarly to the mildly nonlinear 
case, continuous waves appear and their amplitude increases in time (top left). In this 
case the beating does not appear, and the amplitude grows exponentially in time to 
become 10^ times larger than the initial amplitude (top right). The bottom left panel 
for t = 850M shows the wave forms after the “bosenova”. After the amplitude stops 
growing, the gravitational waves with a nontrivial waveform are radiated. Until the 
occurrence of the bosenova, the gravitational wave frequencies are close to twice the 
axion mass, 0 ~ 2/i. But after the bosenova, the gravitational wave frequencies are 
a few times smaller than 2/i. The gravitational wave amplitude approximately shows 
exponential decay in time. After that there appears another phase during which the 
gravitational waves are generated with the frequency ~ 2/i, but keeping a relatively 
larger amplitude compared to the initial amplitude (compare the top left and bottom 
right panels). 

Figure 28 shows the value of a™ for (h, m) = (0,2) (solid line) and (h, m) = (0,0) 
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Figure 27. Snapshots of the real part (solid line, red online) and imaginary part (long 
dashed line, green online) of spectral component Q!q ' of the Teukolsky variable for a 
gravitational perturbation in the logarithmic scale at time t = 400M (top left panel), 
t = 650M (top right panel), t = 850M (bottom left panel) and t = 1200M (bottom 
right panel) in the strongly nonlinear case. 



r* = 200M 



tIM 

Figure 28. Behaviour of real part of the spectral component a® (solid line, red 
online) and (dashed line, green online) of the Teukolsky variable for a gravitational 
perturbation as functions of time observed at the point r* = 200M in the strongly 
nonlinear case. 


(long dashed line) as functions of time observed at the hxed position r* = 200M. We 

( 2 ) 

can conhrm that ^ shows the exponential growth and has quite nontrivial wave form 
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during the bosenova. After that, it decays with frequencies few times smaller than 
2fi. The highly dynamical process of the bosenova generates the m = 0 mode as well. 
Generation of the m = 0 mode reflects the rapid change in the radial distribution of the 
axion cloud, and it shows frequencies different by a factor from 2/x. 


t = 820M 



rJM 


Figure 29. A snapshot of the real part of the spectral components of the 
Teukolsky variable for a gravitational perturbation for fixed rh = 2 and h = 0 (solid 
line, red online), 1 (long-dashed line, green online), 2 (short-dashed line, blue online), 
and 3 (dotted line, purple online) at time t = 820M in the strongly nonlinear case. 


Figure 29 shows a snapshot for t = 820M, but for with m = 2 and n = 0 


(solid line), 1 (long-dashed line), 2 (short-dashed line), and 3 (dotted line). This is the 
time when the bosenova phenomena is about to cease. Similarly to the mildly nonlinear 
case, the higher n modes are also excited and show similar patterns to the h = 0 mode. 

To summarize, in the highly nonlinear case, the gravitational waves show 
quite interesting behaviour. In particular, the gravitational wave amplitude grows 
exponentially until it becomes 10“^ times the initial amplitude. Further, during bosenova, 
gravitational waves with a nontrivial wave forms are radiated with frequencies smaller 
than 2fi. Also, gravitational waves in higher i modes are excited with patterns similar 
to the lower i modes. 


5.3. Discussion on this section 

In this section, by numerical simulations, we have studied the emission of gravitational 
waves sourced by the scalar field in the three cases: the case without self-interaction, 
a mildly nonlinear case, and a strongly nonlinear case. From the results of these 
simulations, the gravitational wave emission from a realistic axion cloud around a 
rotating black hole would be expected as follows. When the scalar field amplitude 
is small, (p 1, continuous gravitational waves are emitted, and as the field amplitude 
grows large, beating appears in them. The bosenova happens at some critical amplitude 
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of scalar field, and associated with it, the amplitude of gravitational waves grows 
exponentially, and then, burst-type gravitational waves with nontrivial waveforms are 
emitted. Under the influence of scalar field self-interaction, gravitational wave forms 
become complicated and may not be modeled by a simple formula. This would make 
it difficult to search gravitational wave signals in the data analysis. But the good news 
is that the self-interaction enhances the amplitude of gravitational waves, even by a 
factor of 10^ during the bosenova, which would be helpful to detect the signals. Since 
bosenovae happen intermittently in a long time scale, burst-type gravitational waves 
would be emitted like a geyser from the axion cloud around a black hole. 

The growth in the amplitude by a factor of 10^ may seem too large. Naively, the 
increase in the gravitational wave amplitude is expected to be a factor of ~ 10^, because 
the field amplitude increases by a factor of 10, and therefore, the increase in the source 
term of the energy-momentum tensor in the Teukolsky equation is ~ 10^. However, in 
the Klein-Gordon case, the source term is 0(10"^) while the emitted gravitational waves 
have amplitude of 0(10“^). This means that the scalar field quasibound states has a 
conhguration to cause a significant phase cancellation. In the strongly nonlinear case, 
the field configuration changes due to the nonlinear self-interaction, and this makes the 
effect of the phase cancellation weaker. This is the interpretation for the significant 
increase in the gravitational wave amplitude during the bosenova. 

In our approach, we calculated the Teukolsky variable for s = —2. For the 
choice s = —2, _ 2 T is related to i/'o of the Newman-Penrose variables. While i/jq is 
very convenient to calculate the ingoing energy flux to the black hole, it is not directly 
related the energy flux formula towards infinity. In order to evaluate the outgoing flux, 
the conversion from 'i/jq to is necessary (see a similar discussion on conversion from 
‘ip 4 to '00 for incoming waves from inhnity for |1H])- Since the relation between the 
amplitudes of -00 and 04 depends on the angular frequency a), we have to calculate 
Fourier transform of emitted gravitational waves, and reconstruct 04. This procedure 
is necessary to predict the observed gravitational wave forms and is under progress. 

Although we do not know the exact waveform, a rough estimate on the detectability 
of the gravitational wave burst emitted during the bosenova can be given as follows. 
Let us consider the case where the black hole candidate, Cygnus X-1, wears a cloud of 
axion scalar field with mass 2.4 x 10“^^ eV, which corresponds to M/i 0.30. Up to 
the bosenova, the gravitational waves show exponential growth with keeping u ^ 2/j, 
that corresponds to ~ 1200 Hz, and there is a period At lOOM during which the 
gravitational amplitude is larger than the Klein-Gordon case by a factor of 10^ or more. 
On the other hand, in our previous paper mi we calculated the amplitude of continuous 
gravitational waves from Gygnus X-1 in the case that it wears the quasibound state of 
a Klein-Gordon field, using the black hole parameters determined by observations (see, 
e.g., US]). Then, the gravitational wave amplitude from the bosenova is obtained just 
by multiplying 10^ to our previous estimate for the Klein-Gordon case. The result is 



(5.12) 
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after substituting the mass M, the distance d, and axion mass /x. Then, multiplying 
the strain amplitude becomes h^ss ~ 10 “^°/a/Hz. This is about 10 times larger 
than the observation threshold of LIGO if the GUT scale decay constant is adopted. 
This shows the possible detectability of gravitational waves from the bosenova. 

Here, we note that the reliability of the estimate here based on the Schwarzschild 
case may be questioned, because there is difference between the bosenovae in the 
Schwarzschild case and the Kerr case. In the Kerr case, the infall of positive energy 
happens as a result of excitation of m = —1 mode, and it is a slower process compared to 
the Schwarzschild case. This may affect the amplitude of gravitational wave burst during 
the bosenova. Therefore, it is necessary to extend our code to the Kerr background case 
in order to make our calculations more realistic and to enable our simulation results 
applicable to observational data analysis. 

6. Summary and discussion 

In this paper, we have examined the three features of time evolution of an axion cloud 
around a rotating black hole by numerical calculations/simulations. In the superradiant 
phase, we calculated the growth rate by the superradiant instability without self¬ 
interaction (i.e., a the massive Klein-Gordon held) including the overtone modes as 
well as the fundamental modes (section]^. Our results indicate that for £ = m = 3, 
the overtone modes can have larger growth rates compared to the fundamental mode 
when the black hole is rapidly rotating. In section we studied the phase where 
the nonlinear self-interaction becomes important by simulating the behaviour of the 
scalar held numerically. The bosenova collapse happens for i = m = 1 mode, and 
it becomes more violent as M/x is decreased. For i = m = 2 mode, our simulation 
results suggest that the bosenova collapse would not happen, but a steady outhow will 
be formed instead. If we consider superposition of two modes, further rich phenomena 
were observed. In section]^ we showed our preliminary simulations of gravitational wave 
emission from an axion cloud in the presence of axion nonlinear self-interaction around 
a Schwarzschild black hole. Although the extension to the Kerr case is necessary, we 
have obtained a lot of indications for these phenomena. In particular, gravitational wave 
burst is emitted when the scalar held collapse to a black hole. More detailed summary 
and discussion are given at the end of each section. 

Up to now, we took attention to the phenomena that happens in the initial 
phase (superradiant instability) or in the time scale of lO^M after that (bosenovae 
and gravitational wave emission). In what follows, we consider the evolution of an 
axion cloud and a black hole in a longer time scale. In a long time scale, the black 
hole parameters M and J change in time due to the backreaction of energy/angular 
momentum extraction by the axion cloud. The evolution of the black hole parameters 
was discussed in the original axiverse papers Dca. and the more detailed adiabatic 
evolution was calculated in [20], taking account of the supply of energy/angular 
momentum by the accretion disk, the energy/angular momentum extraction by the 
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scalar cloud, and the loss of the energy/angular momentum of the axion cloud by the 
gravitational wave emission. When the superradiant instability becomes important as 
the black hole is spun up by the accretion disk, the scalar cloud rapidly grows to decrease 
the black hole spin parameter a* to the marginal value satisfying uj = mfln- After that, 
the spin parameter gradually grows due to the effect of the accretion disk. If we switch 
on the self-interaction effect, the nonlinear effect becomes signihcant when the axion 
held becomes order of the decay constant, $ ~ /a, and the axion cloud energy at this 
time is E/M ~ 10^(/a/Mpi)^. For /a with the GUT scale (1/10 of the GUT scale), 
the corresponding energy is ~ 0.1% (~ 0.001%) of the black hole mass. Therefore, 
the evolution is expected to highly depend on the value of /a, and the scenario would 
become different from the ones in PElEo]. 

One important additional quantity in the presence of the nonlinear self-interaction is 
the energy/angular momentum loss rates of the axion cloud by the bosenova {^ = m = 1 
case) or the axion outhow {i = m = 2 case). Since the direct approach is very difficult 
because the nonlinear self-interaction makes the system highly dynamical, one has to 
consider an effective approximation. If we take an average over the timescale much longer 
than the bosenova occurrence, the energy/angular momentum of the gravitationally 
bounded axion held would be approximately constant, whose values depend only on a*, 
M/i, and /a. The decrease in the energy/angular momentum of the black hole due to 
the superradiant instability in short time scales will be erased out after the averaging 
because of the recycling processes by the bosenova. Then, the system is regarded as a 
converter of the black hole energy/angular momentum into scattered axion helds and 
gravitational waves. Due to this conversion, the black hole decreases the spin parameter 
a* to the marginal value. Because its time scale becomes longer as /a is decreased, we 
may be able to obtain constraints on /a from the fact that some solar-mass black holes 
are rapidly rotating like Gygnus X-1. This would give a method for constraining the 
axion model parameter different from the one of HZ] that uses the direct observation of 
gravitational waves. 

In order to carry out this program, we have to make our gravitational wave code 
applicable to the Kerr case and discuss how to include the backreaction of gravitational 
wave emission to the axion held dynamics. Then, by long term simulations, we will 
be able to determine the averaged conversion rate of the extracted energy/angular 
momentum to the scattered axion helds and gravitational waves. We are planning 
to challenge these problems as the next step. 
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Appendix A. Numerical method: Scalar field 


In this Appendix, we give supplementary informations for the numerical method for the 
scalar field simulations with the pseudospectral approach. 

Appendix A. 1. Equation 

The explicit form for the sine-Gordon equation for (p in the Kerr background spacetime 


IS 


[(r^ + — Aa^ sin^ — AMarip^^ 

A / /I \ 2 k — sin^ 6 2 v^ a • 

+A [q) 0 e + cot H-- (^,^0 - /i SA sm (p. 

sm 6 


(A.l) 


Substituting (4.4) and (4.7), we have 


[(r= + a^f - Ao=(l - y)] /<”> = (r= + a")V,S. + 2rA/;”> 

-(4imMar)/l’") + A(1 - !,2)/<™) - 2(|m| + ^A^/f 

+ [m^a^ — A(m^ + |m|)] — /i^EAs*”*^ (A. 2 ) 


after equating the coefficients of We further substitute (4.6) and (4.8). Equating 
the coefficients of y"', the equation for is found to be 

[(r^ + a^f - = (r^ + + 2rAa^y^l^ 

— (AimMar) + {n + 2){n + l)Aa [^"^2 

+ — (n^ + n + + |m| + 2|m|n) A] — /i^A . (A.3) 

Note that oci the left-hand side must be interpreted as zero for n = 0 and 1. 


Appendix A.2. Nonlinear term 

Here, we discuss how to compute the spectral components for the nonlinear term 
sin(y 9 . This can be done with the Taylor expansion of sin(y 9 . 


sm if = 


E 


(- 1 )^ 

( 2 / + 1 )! 




21+1 


(A.4) 

I Z/, -H M ! 

Z=0,l,2,... 

In the numerical computation, we introduce a cut-off parameter /^ut for 1. The fortunate 
thing is that the convergence radius for this Taylor expansion is infinity. This means 
that by increasing the value of /cut, we can approximate sin(y 9 in a larger range of (p. 


Figure [AT] compares the values of sin a: and its Taylor series up to / = /cut for various 
/cut values. The Taylor series gives a good approximation from zero to a certain value of 
X, and the approximation rapidly becomes bad beyond that point. The range of good 
approximation becomes larger as lent is increased, and the difference is almost invisible 
for /cut = 9 in the range 0 < a; < 7. 

The spectral component of the Taylor series of sin(p can be calculated if we can 
compute each spectral component of For this reason, let us derive the general 
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Figure Al. Comparison between the value of sin a; (solid line) and its Taylor series 
up to the term I — lent for lent = 0, ...,9 (dashed lines). 


relation between the spectral components of two functions and those of their product. 
Let us consider /, g, and h, with 


f = 

E 

sinH 0/(-)(t,r,0)e'™‘^, 

j{-m) 

_ j{m)* 

(1.5a) 


m=0,±l,±2,.. 





9 = 

E 

sinH 0^M(t,r,0)e™'^, 

gi-m) ^ 

^gim).^ 

(1.55) 


m=0,±l,±2,.- 





h = 

E 

sinH eh’^^\t,r,e)e^'^^, 

h(-^) 

_ j^bn)* ^ 

(1.5c) 


m=0,±l,±2,... 


li h = fg holds, we can derive the relation 

m 

hM = E j{k)g(m-k) _ y2'^k (^j{k+m) y(k)* g(k+m) j!(k)*^ ^ 

k=0 k=l,2,... 


Next, let US consider the expansion with respect to y, 

g^^^= and 

n=0,l,... n=0,l,... n=0,l,... 

Then, we hnd the relation 

m n n j 

k=0 j=0 k=l,2,... j=0 i=0 

where we dehned 

EG/2 (i = 0,2,4,...,2fc). 

0 (I = 1,3,5,...), 

0 (i = 2/c + 2,2/c + 4,...). 


(1.7) 


( 1 . 8 ) 


(1.9) 


Using this relation, we can compute iteratively, and thus, the spectral components 

(m) r • 

an oi sm(p. 
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Appendix B. Numerical method: Teukolsky variables 

In this Appendix, we give supplementary informations for the numerical method for the 
gravitational wave simulations. 


Appendix B.l. Equation 


We solve the Teukolsky equation for s = —2 by the pseudospectral decomposition in the 
angular coordinates in a similar manner to the scalar case. The equation for the 

spectral components of is presented in (5.7). We further decompose the coordinate 
y with (5.10). By equating the coefficients of y”, we have 

[(r^ + 0^)2 - + {a^A)a^E\ = (r^ + 

+ 2 (r 2 + a 2 ) 


^ A rA 

A r -1-5^ 

r 

(m) 


~ 4 [M{P — a^) — Ar + imMar] a. 


(m) 

h 


+ (4ia)Adf™; + (^ + 2)(fi + l)Aa‘'^\ - (n + 1) (|m - 2| - |m + 2|) 


(m) 


+ < — 4ima(r — M) + A ( 2-A^^ H-— 


2 A\ 


- 


r2 I 


A 


— {2n^ — 2n + m2 + |m2 — 4| + 3|m — 2| + 3|m + 2|) 


( 2 . 1 ) 


where we have assumed that the source term can be expressed in the form of the series 


dvr 


|m-2| + |m + 2| 

2 2 Lr 


_ y'^\m,-2\/2(^l _|_ ^^|m+2|/2 


^ t 




( 2 . 2 ) 


n=0,l,2,... 


Note that left-hand side must be interpreted as zero for n = 0 and 1. 


Appendix B.2. Regularity of the source term 


In the Schwarzschild case, the source term is reduced to a relatively simple formula 
that guarantees the regularity at the two poles, y = ±1. The explicit form of 



introduce the two differential operators: 


= dt- dr,, 
p( 2 ) ^ (p(i))2 _ 2 (r 

Using these operators, the following formula is derived: 
(4/r2)T^™^’”^2 

(1 _ ^)|mi+m2-2|/2('X + |mi+m2+2|/2 “ ^^2 l"i2|)(l+2/) ^ (1 

+ (m2 - |mi|)(l + y)^il - y)^ fim,) 


(2.3a) 

(2.36) 


y^^Ar- y(™2)p(2) j(mi) 
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p(+) p(-) / A 2 \ 

+2mim2(l + y)^il - y)^ f fim^) ^ \ 

p(+) p(-) / 

+(1 + y)^{l-y)^ + / 7 i)p( 2 )j(m 2 ) _ 2 i)(i)/(-i)p(i)/(- 2 ) + 2 

p(+) p(-) / A2 

+2m2il + y)^ il - y)^ _/72)p{2) j(mi) ^ p(l) j(rni)p(l)^(ma) _ 

\ 5 i/ ,^0 ty 

p(+) p(-) / \ 

+ 2 mi(l +?/)^(l - 7 /)^ _/ 7 i)p{ 2 )y(m 2 ) ^p(l)y(rn 2 )p(l)y.(mi) _ f{mr) ^{^2) ^ ^ 2 . 4 ) 

\ >i/ 5i/ ly^o iy j 


where 

= |mi| + \m 2 \ — 2 + 46'(±m2) — |mi + m 2 ± 2| (2.5a) 

= l"^i| + l"^ 2 | — 2 + 46*(±mi) — |mi + m 2 ± 2| (2.56) 

= l"^i| + l"^ 2 | — 2 + 26*(±mi) + 26*(±m2) — |mi + m 2 ± 2| (2.5c) 

= |mi| + |m2| + 2 - |mi + m 2 ± 2| (2.5d) 

= l"^i| + l"^ 2 | + 26*(±m2) — |mi + m 2 ± 2| (2.5e) 

= |mi| + |m 2 | + 26'(±mi) - |mi + m 2 ± 2| (2.5/) 


Here, all non-negative for arbitrary integers mi and m 2 , and hence, the regnlarity 

of the sonrce term at the two poles is guaranteed. Substituting /("^d = 

/(™' 2 ) = we can derive the expression for for each value of (mi, m 2 ). 


! \ 
_j(mi)y(m2) I 


^ ,2/ ,2/ 
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